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ABSTRACT 

We present a detailed comparison between the well-known smoothed particle hydro- 
dynamics (SPH) code GADGET and the new moving-mesh code AREPO on a number 
of hydrodynamical test problems. Through a variety of numerical experiments with 
increasing complexity we establish a clear link between simple test problems with 
known analytic solutions and systematic numerical effects seen in cosmological sim- 
ulations of galaxy formation. Our tests demonstrate deficiencies of the SPH method 
in several sectors. These accuracy problems not only manifest themselves in idealized 
hydrodynamical tests, but also propagate to more realistic simulation setups of galaxy 
formation, ultimately affecting local and global gas p roperties in the full cosmo logi- 
cal framework, as highlighted in companion papers by Vogelsberger et al. (2011) and 
Keres et al. (2011 ). We find that an inadequate treatment of fiuid instabilities in GAD- 



GET suppresses entropy generation by mixing, underestimates vorticity generation in 
curved shocks and prevents efficient gas stripping from inf ailing substructures. More- 
over, in idealized tests of inside-out disk formation, the convergence rate of gas disk 
sizes is much slower in GADGET due to spurious angular momentum transport. In 
simulations where we follow the interaction between a forming central disk and orbit- 
ing substructures in a massive halo, the final disk morphology is strikingly different 
in the two codes. In AREPO, gas from infalling substructures is readily depleted and 
incorporated into the host halo atmosphere, facilitating the formation of an extended 
central disk. Conversely, gaseous sub-clumps are more coherent in GADGET simula- 
tions, morphologically transforming the central disk as they impact it. The numerical 
artifacts of the SPH solver are particularly severe for poorly resolved fiows, and thus 
inevitably affect cosmological simulations due to their inherently hierarchical nature. 
Taken together, our numerical experiments clearly demonstrate that AREPO delivers 
a physically more reliable solution. 
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1 INTRODUCTION 

Numerical simulations have become an indispensable tool for 
studying astrophysical phenomena. Over the last decade, the 
numerical accuracy and fidelity of simulation methods have 
undergone drastic improvements, both in terms of resolvable 
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dynamic range and of the complexity of the physical pro- 
cesses that are now routinely incorporated into the codes. 
This remarkable progress in numerical techniques coupled 
with the ever increasing power of high performance com- 
puting platforms has led to major advances in a number of 
topics, such as star formation ( Klessen et al.||2009|, accre- 
tion disk dynamics and associated jet phenomena (|Gammie 



et al. 2003 De Villiers & Hawley 2003), supernova explo- 
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sions ([Janka et al.| |2007 ), black hole coalescence fPretorius' 
20071 , and cosmic structure formation (Springel et al. 2005). 



discrepancy as found by Frenk et al. (1999), which they at- 



Clearly, an in-depth understanding of many astrophys- 
ical processes, particularly those that are inherently non- 
linear, relies crucially on the accuracy and realism of the 
numerical approach. While the latter is largely determined 
by the degree of adequateness, comprehensiveness and con- 
sistency of the physical model assumed, the former depends 
sensitively on the discretization scheme adopted for the 
equations, which with increasing resolution should converge 
to the correct solution. Nonetheless, even for some elemen- 
tary physical problems, with known analytic solutions, sim- 
ulation methods can sometimes produce poorly converged 
results, or even converge to a solution which is, however, 
different from the expected result | Springel||2010b ). 

In the context of hydrodynamical cosmological simula- 
tions, an important example is given by the seminal work of 



Frenk et al. ( |1999J (the Santa Barbara Comparison Project). 
This study performed a detailed comparison between twelve 
different simulation codes that tracked the non-radiative 
evolution of a forming galaxy cluster in a cold dark mat- 
ter (CDM) cosmology. The initial conditions were gener- 
ated independently by each group from a provided linear 
theory density or displacement field. Very different numeri- 
cal methods, which can be broadly classified into smoothed 
particle hydrodynamics (SPH) and mesh-based techniques, 
employing also different gravity solvers and different effec- 
tive resolutions, were then used by the groups to follow the 
formation and evolution of the target object. Reassuringly, 
the Santa Barbara Comparison Project has shown that the 
global properties of the simulated object, both in terms of 
dark matter and gas distribution, were in reasonable agree- 
ment among the codes. However, detailed gas properties of 
the Santa Barbara cluster, notably in the central region, 
exhibited a much poorer level of consistency. In particular. 



Frenk et al. ( jl999J pointed out that there is a systematic dis- 
crepancy in the central entropy profiles predicted by SPH 
and mesh-based codes, with the former producing power- 
law entropy profiles all the way to the centre and the latter 
yielding some form of entropy core. 

Even though progress in numerical techniques has led to 
substantial improvements in current simulation codes com- 



pared to those considered in Frenk et al. (1999), the sys- 



tematic difference in the predicted central cluster entropy 
still persists today, and its origin has not been fully under- 
stood thus far. Agertz et al.' ("2007) have performed a set 
of numerical experiments with SPH and mesh codes, where 
they have simulated the evolution of a cold, dense blob in a 
hot windtunnel. By comparing the outcomes from different 
methods (and at various resolutions) against the character- 
istic blob disruption timescale expected analytically, they 
conclude that in SPH codes the development of fluid in- 
stabilities can be numerically hampered. The suppression 
of dynamical fluid instabilities, such as Kelvin-Helmholtz, 
Ray leigh- Taylor, and Richtmyer-Meshkov instabilities, leads 
to less efficient mixing of fluid elements with different specific 
entropy, and hence also inhibits entropy generation through 
mixing in the simulated system. In simulations of colliding 



isolated galaxy clusters, Mitchell et al. (2009) have shown 



that the central entropy profiles obtained with the SPH code 
GADGET and the mesh code FLASH show a similar level of 



tribute to the different levels of mixing. 

There have been a number of attempts to improve the 
description of mixing in SPH by modifying the standard 
set of discretized equations (see e.g. 



Price 



2008 



Wadsley 



et al.|2008 HeB &: SpringelpOlOf or by incorporating a Rie- 



mann solver in place of the artificial viscosity (e.g. Inutsuka| 



2002| |Cha Whitworth|[2003l [Murante et al.||2011t . For 



example. Price (2008) has shown that the introduction of 



an artificial thermal conductivity term can improve the be- 
havior of fluid elements at contact discontinuities, which in 
turn leads to better developed Kelvin-Helmholtz instabili- 
ties. More generally, Wadsley et al.| (|2008) suggested that a 
physical modeling of heat diffusion is necessary when simu- 
lating high Reynolds number flows (both for SPH and mesh- 
based methods), and that by applying such an approach a 
flat entropy core is likely a more correct solution for non- 
radiative galaxy cluster simulations. 

Differences in the hydrodynamical solver between SPH 
and mesh-based codes are possibly not the only reason for 
the systematically different central entropy profiles in the 
Santa Barbara Comparison Project. As mentioned above, 
the groups involved in this study did not use identical initial 
conditions, and did not perform their simulations at equal 
numerical resolutions and with identical gravity solvers, 
which opens up the possibility of additional sources for dis- 
crepancies. A more recent code comparison study by |Heit-] 



mann et al. (2008) evolved uniform, dark matter only cos- 



mological boxes in a ACDM universe with 10 different codes, 
starting from the same initial conditions. They showed that 
for large systems there is a reassuring agreement in the halo 
mass functions between the codes, as well as in the internal 
structures of haloes in the outer regions. However, the study 
by Heitmann et aL] (|2008) also revealed significant discrepan- 
cies for small halos, demonstrating that the typical root grid 
resolution commonly adopted in adaptive mesh refinement 
(AMR) codes for simulations of cosmic structure formation 
is overly coarse and leads to a suppression of low mass halo 
formation (see also O'S hea et al.|2 005). This emphasizes the 
need to simultaneously strive for high accuracy both in the 
gravity solver and in the hydrodynamics. Additionally to the 
high accuracy of the gravity and hydro solver, it is also of 
prime importance to adopt sufficiently high mass and spatial 
resolution for the studied problem at hand. For example, if 
besides pure hydrodynamics other physical processes, such 
as star formation and associated feedback, are modelled in 
cosmological simulations, it is necessary to resolve all star- 
forming haloes with a sufficiently large number of resolution 



elements to reach convergent results (| Springel Hernquist 



2003b ). 

In the present work we perform a detailed comparison 
study between the widely used SPH code GADGET and the 
novel moving-mesh code AREPO on a number of hydrody- 
namical test problems with increasing levels of complexity. 
We have adopted a combination of existing test problems 
and newly devised numerical experiment^ in order to pro- 
vide a clear link between the results of our test problems and 



High-resolution images and movies of various numeri- 
cal experiments are available for download at the website 



|http://www. cfa.harvard.edu/itc/research/movingmeshcosmology 
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those of full cosmological simulations, which are discussed 
in detail in our companion papers ( Vogelsberger et aL|2011| 
hereafter Paper I) and ( [Keres et al.|2011 ^ hereafter Paper II). 
Our work thus contributes to the understanding of the signif- 
icant differences in baryon properties found in cosmological 
simulations between mesh-based and SPH techniques, both 
at a global level (see Paper I) and at the scale of individual 
galaxies (see Paper II). 

A unique advantage of our comparison study lies in the 
fact that simulations with GADGET and AREPO can be 
started from identical initial conditions and that both codes 
employ the same gravity solver. Note that in case of AREPO 
gravitational softenings for the gas can be kept fixed as is the 
case of GADGET or can be computed adaptively determined 
by the cell size. In the numerical experiments with gas self- 
gravity we explore both fixed and adaptive gas softening, 
where in the case of adaptive softenings a floor equal to the 
GADGET gas softenings is set. We find that these different 
choices of gravitational softenings in AREPO do not affect 
our findings. This allow us to isolate cleanly how the dif- 
ferent hydro solvers used by GADGET and AREPO impact 
gas properties. We first consider elementary hydrodynami- 
cal numerical tests such as a strong ID Sod shock tube test, 
a 2D implosion test (which is widely used for benchmark- 
ing mesh codes, but has rarely been considered in SPH), 
and the "blob" test ( |Agertz et al.|2007 ). We then perform a 
number of isolated or merging halo simulations in the non- 
radiative regime, aimed at understanding how shocks and 
fluid instabilities affect their gaseous atmospheres. Finally, 
we study the differences between GADGET and AREPO in 
radiative simulations, where we follow inside-out disk forma- 
tion and interactions between the central disks and orbiting 
substructures. 

This paper is organized as follows. In Section [2] we pro- 
vide a brief overview of the numerical codes used and the 
types of simulations performed. Section[3]represents the core 
of the paper where all of our numerical experiments are dis- 
cussed. Finally, we summarize our findings in Section [4] 



2 METHODOLOGY 
2.1 Numerical codes 

2AA GADGET 

GADGET is a massively parallel TreePM-SPH code widely 
used in numerical astrophysics. In this study we adopt the 
latest GADGET-3 version. A detailed description of an ear- 



lier version can be found in Springel ( 2005 ). GADGET is fully 
adaptive in time and space, and in its entropy formulation 
for SPH (Springel & Hernquist 2002) manifestly conserves 
both energy and entropy in the absence of artificial viscosity. 
Gravitational forces are computed with an octtree method 
( Barnes Hut [1986 Hernquist|1987 ). To speed up the com- 
putation, long-range forces can be optionally evaluated with 
a PM method, with the tree being restricted to short-range 
forces only. 

In all our tests, we adopt as standard value for the artifi- 
cial viscosity strength a = 1.0, and we use 64 neighbours for 
kernel interpolation in three dimensional simulation, unless 
we specifically vary these parameters to assess their effect. 



2.1.2 Other SPH implementations 

As mentioned in the introduction, there have been a num- 
ber of recent proposals to improve the standard SPH imple- 
mentation in various ways, for example by invoking a time- 



dependent artificial viscosity ( ^Morris Monaghan 1997 



Dolag et al. 12005 ), a modified density est imate (e.g. [Ritchie 
Thomas] [2001] |HeB Springel] [^hlQl ]Saitoh M^km^ 
2012|, a decoupling of the hot and cold neighbours in mul- 
tiphase flows ( ]Marri k Whit e 2003 j Oka moto et al.|]2003| ), 
an artificial therm al conductivity (,Pricejj2008 ), a modified 
equation of motion ( |HeB Springel 2010 ]Abel 2011 ), an ex- 
plicit modelling of mixing ("Wadsley et al.j2008), an enlarged 



neighbour number combined with a different kernel shape 
( ]Read et al.|2010| ), or a replacement of the artificial viscosity 
by a Riemann solver (e.g. Inutsuka|2002 Cha & Whitworth 



2003| ]Murante et~aLl]2011| ). While some of these modifica- 
tions of the standard SPH formalism deliver more accurate 
results in targeted numerical experiments, no consensus has 
yet emerged whether any of these approaches (or a combina- 
tion thereof) is sufficiently robust for cosmological applica- 
tions and leads to universally more accurate results in galaxy 
formation simulations. Therefore, in this work we focus on 
the traditional SPH implementation rather than on the var- 
ious possible modifications proposed recently. We note that 
this standard formulation of SPH has also been typically 
employed in state-of-the art cosmological SPH calculations 
(e.g. Grain et al.|2009l ]Di Matteo et al. 2012). Furthermore, 
as discussed in Paper I, there are other, generic issues with 
SPH that have not been resolved by any of the above mod- 
ifications. These issues ultimately mean that SPH does not 
currently have a formal convergence condition, which also 
complicates rigorous evaluations of variants of the standard 
SPH algorithm. 



2.1.3 AREPO 



AREPO (Springel 2010a) is a new massively parallel sim- 



ulation code, which uses the same gravity solver as GAD- 
GET (augmented with the possibility of adaptive gravita- 
tional softenings for the gas), but employs a completely 
different method for the evolution of the fluid. It adopts 
a second-order accurate finite volume technique, where the 
solution of the Euler equations is computed by an unsplit 
Godunov method equipped with an exact Riemann solver. 
Throughout we use the default choice of the slope limiter 
in AREPO which prevents the linear reconstruction to over- 
or undershoot the maximum/minimum values of neighbour- 
ing cells, as described in detail in "Springel' f2010a). Unlike 
standard finite- volume codes used in numerical astrophysics, 
AREPO solves the equations on an unstructured Voronoi 
mesh, which is allowed to freely move with the fluid. The 
resulting quasi-Lagrangian nature of AREPO automatically 
guarantees spatial adaptivity and greatly reduces numerical 
diffusivity even in the presence of large bulk flows. Com- 
pared to standard Eulerian mesh codes, AREPO has the ad- 
vantage of being fully Galilean invariant (as is GADGET as 
well), it is less prone to advection errors and over-mixing, 
and preserves contact discontinuities better. 
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In this study, we employ a mesh regularization metho(|^ 
in AREPO by default, based on a Lloyd algorithm (for de- 
tails see Springel |2010a| , which ensures that the geometric 
centre of each Voronoi cell is sufficiently close to the cell's 
mesh-generating point to ensure good accuracy of the spa- 
tial reconstruction. In cosmological simulations (see Paper 
I), an alternative regularization criterion has proven to be 
advantageous, based on the maximum opening angle under 
which a cell face is seen from the mesh- generating point. 
We have checked for a number of test runs presented in this 
study that this alternative regularization method does not 
affect any of the results described here. 

For most of the simulations presented here we do not 
use the possibility of mesh refinement and de-refinement op- 



erations (see Springel 2010a), except in our numerical ex- 



periments with star formation presented in Section |3.4.3| 
where we employ it for verification of our findings. Mesh 
de-/refinements are used to constrain the mass of cells to a 
small range around a target value (equal to the gas parti- 
cle mass in the matching GADGET run). This restricts the 
spectrum of star particle masses which are generated from 
gaseous cells, and thus assures that N-body heating effects 
are minimized. Also, as our default choice we use the en- 
ergy formulation of AREPO. We verified for each numerical 
experiment that there is no significant spurious transfer of 
kinetic into thermal energy. 



2.2 Types of simulations 

2.2.1 Physical processes 

We perform both radiative and non-radiative hydrodynam- 
ical simulations, where in the former case gas is represented 
by a primordial mixture of hydrogen and helium in an op- 



tically thin limit (Katz et al. 1996). In simulations with 
radiative cooling, we employ a subresolution multi-phase 
model for star formation and associated supernova feedback 



(Springel & Hernquist 2003a). We slightly modify the be- 



haviour of this model for gas elements which are hot but 
already above the density threshold for star formation, by 
allowing them to settle quickly onto the effective equation 
of state: if their newly estimated temperature from radiative 
cooling would fall below the temperature of the multiphase 
medium, we set it equal to the multiphase temperature. This 
change has been introduced to make the sub-grid star for- 
mation module consistent with its current implementation in 
AREPO. We also perform some simulations with the subreso- 
lution star formation model where spawning of star particles 
is intentionally prevented, but the cold, dense gas above the 
density threshold for star formation is still governed by the 
effective equation of state. These simulations are particu- 
larly useful for understanding the development of dynami- 
cal instabilities between cold and hot media, and have the 
advantage of not being prone to fragmentation which might 
affect pure cooling runs. 

Note that we deliberately consider only very simplified 
baryonic physics implementations in all presented numeri- 
cal experiments, as detailed above. This is for two reasons. 



^ Note that in the presence of mesh regularization the Galilean 
invariant nature of the AREPO code is not violated. 



First, we need to make sure that gas cooling and star forma- 
tion processes are treated on an as equal footing as possible 
in order to allow a meaningful comparison of the two codes. 
Thus, we adopt a simple multi-phase model for star forma- 
tion which is largely insensitive to the detailed structure 
of the gas on very small scales where gas joins the inter- 
stellar medium and is described by a relatively stiff effec- 
tive equation of state. Second, our aim is to isolate in an 
as clean manner as possible the differences between simu- 
lated systems with GADGET and AREPO stemming from 
the discretization of the fluid equations alone. This goal 
is given precedence over trying to reproduce observational 
findings through a more sophisticated modelling of addi- 
tional physics. While it is likely that such additional physi- 
cal processes will change the properties of some of our simu- 
lated systems by possibly different degrees in GADGET and 
AREPO, it is of significant interest in its own right to dis- 
entangle numerical uncertainties from uncertainties in the 
physical modelling of star formation and associated feed- 
back processes. 



2.2.2 Gravitational softenings 

For simplicity, many of the numerical experiments presented 
in this study have been evolved without gas self-gravity. For 
simulations where self-gravity of the gas is nevertheless in- 
cluded, we note that GADGET employs constant gravita- 
tional softening for gas particles in the manner of Hernquist] 



& Katz (1989), while AREPO uses either the same constant 



gravitational softening or an adaptive softening determined 
by the cell size with a minimum softening value set to the 
fixed softening of the matching GADGET run. We have veri- 
fied that this does not lead to substantial differences for any 
of the tests presented in this study. 



2.2.3 Effective hydro resolution 

Even though GADGET and AREPO calculations use the 
same gravity solver and can be initiated from identical ini- 
tial conditions, due to the very different nature of the hydro 
solver it is not straightforward to define unambiguously a 
comparison strategy at the same or equivalent hydro reso- 
lution. In this study we choose to keep the number of res- 
olution elements the same in the both codes, i.e. to have 
the same number of SPH particles as Voronoi cells, corre- 
sponding to a comparable mass resolution in the gas. This 
allows us to adopt indeed the same initial conditions and to 
also have similar mass resolution in the stellar component 
in those simulations where star formation is included. Fur- 
thermore, the CPU costs of the two codes are then roughly 
comparable, as discussed in Paper I. 

We note, however, that this choice implies that the "ef- 
fective" spatial resolution of GADGET is lower than that 
of the matching moving-mesh calculation, given that fluid 
properties in SPH are evaluated by kernel averaging over 
a somewhat larger number of neighbours than needed in 
AREPO for its stencil of neighbouring cells, for example for 
gradient estimates. For this reason, we perform all of our 
numerical experiments at a number of different resolutions, 
which also help us to gain some insight into the convergence 
properties of the two codes. Nonetheless, it is important to 
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stress that the convergence properties of the SPH method 
are still not well understood. For example, one would ulti- 
mately require that the number of neighbours be increased 
with increasing total particle number (see e.g. Rasio|[2000 



[Read et aL]|2010| [Robinson Monagh'aii||2011[ and discus- 
sion in Paper I) , but the appropriate scaling of the neighbour 
number with increasing resolution is currently unknown. It 
is common practice, which we adopt as well, to simply al- 
ways keep the number of neighbours fixed when the total 
particle number is varied, even though the discretized rep- 
resentation of the density field is not guaranteed to converge 
to its underlying smooth distribution with this choice. 



2.3 Initial conditions generation 

For all the tests presented in this study (except for the 
"blob" experiment of |Agertz et al.|[2007 , for which we take 
publicly available initial conditions) we generate the initial 
conditions in terms of simple particle/point distributions. 
Depending on the simulated problem at hand, we either pop- 
ulate the simulated domain with particles uniformly spaced 
on a regular lattice or we adopt Poisson sampling of a den- 
sity field. For simulations with AREPO, the positions of par- 
ticles in the initial conditions define mesh generating points. 
All purely hydrodynamical numerical experiments have been 
performed with periodic initial conditions. In the case of sim- 
ulated problems where we follow the evolution of an isolated 
gaseous halo or merging haloes, we select vacuum boundary 
conditions for GADGET. In simulations with AREPO, it is 
however necessary to enclose the simulation in a finite vol- 
ume in order to make the tessellation with a mesh well posed. 
If needed, we thus add a low resolution background grid to 
our initial conditions, which is chosen sufficiently large as 
to encompass the spatial extent of our simulated objects at 
all times. The procedure adopted for this background grid 
generation, which also reduces the Poisson noise in the ini- 
tial mesh geometry, is described in detail in Section 9.4 of 



Springel (2010a) 



3 RESULTS 

3.1 Strong shocks and interacting curved shocks 
in multi-dimensions 

3.1.1 Strong shock in ID 

As an introductory problem, we consider a strong shock with 
a Mach number M = 6.3 in one dimension. The initial con- 
ditions have been constructed as follows: in a computational 
domain of length Lx = 10.0, N = 200 particles (cells) have 
been placed on a regular grid such that for x < 5.0 the pres- 
sure and density are P = 30.0, p = 1.0, while for x > 5.0 
they are P = 0.14 and p = 0.125, respectively. We adopt 
an adiabatic index 7 = 1.4, and assume that the fluid is 
initially at rest. In the test run with GADGET, we adopt a 
standard value of the artificial viscosity equal to a = 1.0, 
and we vary the neighbor number A^ngb by setting it to 5, 7, 
11, or 15, appropriate for the ID nature of the test. 

In Figure [l] we show the gas density, velocity, en- 
tropy (i.e. P/p^) and pressure at time 0.25 for GADGET 
with A^ngb = 5 (left-hand panels), and for AREPO (middle 



panels). Blue symbols give the values of individual parti- 
cles/cells, dashed red lines represent the initial conditions, 
while dotted red lines are the analytic solution to this Rie- 
mann problem. It can be seen that in both GADGET and 
AREPO the post-shock properties of the fluid are captured 
well, but the shock and the contact discontinuity are signifi- 
cantly broader in GADGET which also shows a characteristic 



'pressure-blip' at the contact discontinuity (see also [Springel 
2010b). Contrary to what one may perhaps suppose, the 



post-shock oscillations present in GADGET are not caused 
by inadequate artificial viscosity, but are instead induced by 
an inaccurate treatment of the sharp contact discontinuity 
of the initial conditions. To demonstrate this point, we have 
performed exactly the same shock tube test but this time 
smoothing the initial contact discontinuity over 5 particles 
in density and internal energy with a Hann window func- 
tion, so as to reduce its sharpness. Green cross symbols in 
the third row of Figure [l] (see also the right-hand panels 
where we zoom into the region around the shock) illustrate 
how the gas entropy is affected by this choice of smoothed 
initial conditions. Post-shock oscillations and the "spike" in 
entropy at the contact discontinuity are greatly reduced. 

It is interesting to note that also in the case of AREPO 
the "spike" in the entropy at the contact discontinuity is 
cured by our smoothed initial conditions. This feature is 
absent in the results of standard grid codes, as we have veri- 
fied by running this test with a static mesh option in AREPO 
(and by running it with ENZO), which tends to broaden con- 
tact discontinuities (depending on the advection speed) and 
thus largely wash-out this feature. This can be seen in the 
lower right-hand panel of Figure^ where we show the results 
from the AREPO run with a static mesh with magenta trian- 
gle symbols. The much lower numerical diffusivity of AREPO 
for contact discontinuities preserves the initial start-up fea- 
ture to much higher accuracy, but at the same time this can 



lead to a larger 'wall heating' effect (Rider 2000) than in 
static grid codes, which tend to smooth out at some level 
the initial start-up errors at the contact discontinuity (see 
also the description of the Noh problem in Springe l||2010a ). 

When we increase the number of neighbors from iVngb = 
5 (as shown in Figure [T]) to 15 in the shock tube tests per- 
formed with GADGET (but keeping the total number of par- 
ticles constant), the spatial region over which the contact 
discontinuity and the shock are broadened increases pro- 
gressively, but the pressure jump between the post-shock 
and pre-shock gas outside of the broadened region remain 
the same, yielding consistent Mach numbers. Based on these 
tests we later discuss in more detail the consequences of 
shock broadening in GADGET when simulating the radial 
infall of gas into static dark matter haloes in Section [3. 3. 2| 



3.1.2 Interacting shocks in 2D 

While one-dimensional shock tube tests are an essential ba- 
sic benchmark for hydrodynamical code performance, they 
are far less demanding than multi-dimensional flow where 
complex interactions of non-planar shocks may occur, as is 
the case in realistic structure formation simulations. Such 
tests have however not been examined widely in the SPH lit- 
erature thus far. We hence perform the so-called "implosion 
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Figure 1. One dimensional shock tube problem with Mach number M = 6.3 at time t = 0.25. Left-hand panels: GADGET with standard 
artificial viscosity a = 1.0 and A^ngb = 5. Middle panels: AREPO with moving mesh. Dashed red lines denote initial conditions, while 
dotted red lines represent the analytic solution. In the third row, green cross symbols illustrate gas entropy in a shock tube of the same 
strength but with initially smoothed contact discontinuity adopting a Hann window. Right-hand panels: zoom into the entropy profile 
around the shock - here magenta triangle symbols in the lower panel are for the AREPO run with a static mesh. 



test" in two dimensions (" Hui et al.|l999l p] To set up this test 
problem, we select a computational domain Lx = Ly = 0.3 
with A^x = Ny = 200 particles or cells, respectively (we have 
also run higher resolution versions with A^x = Ny — 400 
and 800). The initial pressure and density are P = 1.0 and 
p = 1.0 for x+y > 0.15, while P = 0.14 and p = 0.125 other- 
wise. The adiabatic index is 7 = 1.4 and the fluid is initially 
at rest. We performed this test with GADGET adopting an 
artificial viscosity of a = 1.0 and A^ngb = 22 smoothing 
neighbours (suitable for this 2D test). While previously this 
test has been considered using reflective boundary condi- 
tions, we here adopt periodic boundary conditions due to 
their more straightforward implementation in SPH codes. 

In Figure [2] we show density maps for simulations per- 
formed with GADGET (left-hand panels) and with AREPO 
(right-hand panels), at times t = 0.1, 0.3, 0.5 and 0.7 



for 



200. The complex, evolving gas density 



See also "Comparison of Several difference schemes 
on ID and 2D Test problems for the Euler equa- 
tions" by Liska, R. and Wendroff, B., on |http:/ /www-| 
itroja.fjfi.cvut.cz/^liska /CompareEuler/compare8/j^ and 



structure is caused by the continuous interaction of shocks 
throughout the computational domain. Initially, due to the 
discontinuity in density and pressure along x -\- y = 0.15, a 
planar mild shock front develops perpendicular to the x = y 
diagonal traveling towards the origin. Given that we have 
adopted periodic boundary conditions, the gas will interact 
on all four sides of the simulated box, and in particular, in- 
teracting shock waves in the lower left-hand corner result 
in the formation of a narrow jet along the x — y diago- 
nal (which is clearly visible only in the AREPO run). As the 
traveling shocks accelerate the fluid in the regions of density 
discontinuity a Richtmyer-Meshkov instability develops, as 
manifested by "mushroom-like" features in the gas distribu- 
tion. 

While the global gas density distribution in Figure |2] 
agrees reasonably well in GADGET and AREPO runs (i.e. 
with respect to the shape of the regions with different den- 
sities and the magnitude of the density differences) , there are 
several significant differences: i) shocks and contact discon- 
tinuities are much sharper in AREPO, in agreement with our 
findings in Section ; 



3.1.1 



http://www.astro.princeton.edu/~jstone/Athena/tests/ 



a) the density distribution in GAD- 
GET appears not only smoothed out, but it is also noisier, 
as evidenced by the graininess of the density maps, which is 
caused by intrinsic noise in multidimensional flows in SPH 
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Figure 2. Implosion test in 2D at times t = 0.1, 0.3, 0.5 and 0.7. Left-hand panels: GADGET (a = 1.0, A^ngb = 22). Right-hand panels: 
AREPO moving-mesh. For each row, the density scale is the same in the left-hand and right-hand panels, covering the following density- 
range: pgas = 1.2 - 0.4. 
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Figure 3. Vorticity maps of the implosion test in 2D at time t = 0.3. Upper panels: GADGET runs at two different resolutions with 
A^x = Ny = 200 (left-hand panel) and A^x = A^y = 800 (right-hand panel). For both resolutions, vorticity maps are computed by finite 
differencing the velocity field. Lower panels: vorticity maps of the AREPO run with A^x = A^y = 200 computed by finite differencing 
the velocity field (left-hand panel), and by calculating vorticity in the code based on a discretized curl-operator directly applied to the 
Voronoi cells (right-hand panel). 



( SpringelpOlOb ); in) most strikingly, perhaps, is the quite 
different appearance of the low density gas in the bottom 
left corner - a narrow, extended, dense jet along the di- 
agonal is largely absent in GADGET (due to the broaden- 
ing of the contact discontinuity) and Richtmyer-Meshkov 
instabilities along discontinuities are suppressed. We have 
verified that higher resolution simulations with GADGET 
(A^x = Ny = 400 and 800) result in sharper shocks and 
contact discontinuities (with a somewhat feeble jet form- 
ing), but they cannot satisfactorily cure the absence of 
well-developed Richtmyer-Meshkov instabilities. This is in 



agreement with previous studies (e.g. Agertz et al. 2007 



Springel 2010b) indicating fundamental limitations of the 
SPH method in its widely used form. We also note that the 
moving mesh in AREPO does an excellent job of maintain- 
ing symmetry across the x = y diagonal. Furthermore, the 
level of numerical diffusion of discontinuities is very low, as 
indicated by the narrowness and length of the jet. 

In addition to the density fields it is also instructive 
to analyze the vorticity distribution in the implosion test. 



as shown in Figure |3] To construct vorticity maps we ei- 
ther i) first compute spatially adaptive velocity field maps 
for each component on a uniform grid (in our case only x 
and y components) and then finite difference them to obtain 
the curl or ii) we compute spatially adaptive vorticity maps, 
where the curl has been evaluated directly in the code. For 
GADGET runs, vorticity maps computed with method i) are 
shown in the upper panels of Figure [3]for = Ny = 200 
and = Ny = 800. Using method ii) vorticity maps look 
essentially the same. In the case of AREPO (shown in the 



lower panels for the same run with A^x 



200) the 



vorticity maps are also very similar when computed with 
method i) or ii). However, regardless of the exact details of 
the vorticity map generation there are marked differences 
in the vorticity field between GADGET and AREPO. For 
the N^ = Ny = 200 run, the vorticity map in GADGET is 
largely featureless and noisy, with artificial suppression of 
vorticity generation at locations where surfaces of constant 
density and of constant pressure are not aligned (i.e. baro- 
clinic source term oc Vp x Vp). Strikingly, even in the case 
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of the high resolution GADGET run with ^ Ny ^ 800 
(upper right-hand panel) vorticity generation in the regions 
with large baroclinic term is largely suppressed compared to 
the AREPO simulation, even though the latter has a much 
smaller number of resolution elements. It is also worth not- 
ing that the overall magnitude of the vorticity is relatively 
high in GADGET with respect to AREPO in the regions 
which should be characterized by having low vorticity. This 
is due to the noise alluded to above, visible both in the den- 
sity and vorticity maps. This noise is inherently present in 
multidimensional SPH simulations. It is caused by inaccu- 
rate pressure gradient estimates and is particularly evident 



in the subsonic regime (for further details see e.g. Springel 
2010b| [Bauer k Springel|[20T2t . The jitter in gas velocities 
caused by the noisy gradient estimates introduces an artifi- 
cial vorticity 'floor' throughout the simulated volume of the 
box. The poor treatment of vorticity in GADGET has also 
direct consequences for the effectiveness of the widely used 
Balsara switch ( Balsara]|1995 ) for artificial viscosity which 
relies on evaluation of velocity divergence and curl. 

Hence, from the results of the "implosion test" we con- 
clude that while GADGET can accurately capture the main 
fluid properties in the case of dynamically interacting shocks 
with complicated geometries, there are systematic biases in 
detailed aspects, related to the development of fluid insta- 
bilities and to the level of fluid mixing at different entropies. 
These unwanted features lead to the suppression of angu- 
lar momentum transport by vortices and to the damping of 
turbulence in the wake of curved shocks (see also Bauer & 
Springel|2012t . 



3.2 Dynamical fluid instabilities 

3.2.1 Blob experiment 



Our findings from the implosion test (Section 3.1.2 ) indicate 
that there are fundamental differences in the treatment of 
fluid instabilities between GADGET and AREPO. To investi- 
gate this issue in detail, we first consider the so-called "blob" 



test as proposed by.Agertz et al. (2007), which has been an- 



alyzed with many different codes by now ( 


Agertz et al.|2007 


HeB k Springel|2010 Murante et al.|2011 


). The idea behind 



this test is to simulate the evolution of a dense cold blob 
immersed in a hot windtunnel, mimicking in a simplified 
way the motion of a satellite galaxy through the intraclus- 
ter medium (see also | Hefi k Springel|[2011| . If the relative 
velocity between the blob and the surrounding medium is su- 
personic, a bow shock will develop in front of the blob. Addi- 
tionally, dynamical instabilities (mostly of Kelvin-Helmholtz 
and Ray leigh- Taylor type) will grow in the subsonic part of 
the flow between the bow shock and the surface of the blob 



(see Agertz et al. 2007 for a detailed description). These 
instabilities will greatly influence the evolution of the blob, 
leading to its eventual disintegration on the characteristic 
Kelvin-Helmholtz timescale set by the initial conditions. 

To simulate this problem, we adopt the initial con- 
ditions used in the original Age rtz et al.| ([2007) paper, 
which are publicly availabl^ The simulated domain con- 
sists of a periodic box with extensions = 2000 kpc. 



Ly = 2000 kpc, Lz = 8000 kpc, and the blob is initially 
placed at Xc = 2/c = = 1000 kpc. The blob has a 
radius Ruoh = 197 kpc, and it is ten times colder and 
denser than the surrounding medium, such that pressure 
equilibrium is ensured. The density and temperature of the 
external medium is Pmedium = 3.13 x lO^^Mokpc"^ and 
Tlnedium = 10^ K, respectively, and it is moving with a 
constant velocity t'medmm = 1000 km s~^ along the z-axis. 



The adiabatic index is set to 7 = 5/3. Following Agertz 



et al. (2007), we define the characteristic Kelvin-Helmholtz 
timescale as txH = 1-6 tcr, where tcr is the blob crushing 
time defined as tcr = 2i?biob(pbiob/pmedmm)^^^/'i^medmm. For 
these initial conditions, the characteristic Kelvin-Helmholtz 
timescale is txH ^ 1.98 Gyr. 

We have performed the blob test at three different res- 
olutions with GADGET and AREPO. The resolutions used 
are: 32 x 32 x 128, 64 x 64 x 256, and 128 x 128 x 512. The 
two higher resolutions correspond exactly to the resolution 
of the simulations used in Agertz et al. ( 2007 ) , while we con- 



structed the initial conditions for the lowest resolution run 
by sub-sampling. 

In Figure |4] we show projected surface density maps of a 
thin slice {Ay — 100 kpc) centred around the blob position. 
Left-hand panels illustrate GADGET runs at times t — txH, 
t = 2tKH and t = 3tKH, while the right-hand panels are 
for the simulations with AREPO in the moving mesh mode. 
The position and the shape of the bow shock are rather sim- 
ilar between GADGET and AREPO runs, especially at early 
times, but the shock is broader and less crisp in GADGET. 
In the run with GADGET, the blob acquires a cap-like ap- 
pearance caused by the internal shock which compresses it, 
and it undergoes continuous ablation due to the low pres- 
sure region which forms in the wake of the blob ([Agertz 



http://www.astrosim.net/ 



et al. [2007). While initially the blob evolution is similar in 
AREPO, after t — txH well developed Kelvin-Helmholtz and 
Ray leigh- Taylor instabilities lead to an efficient shredding of 
the blob, mixing it with the external medium. In Figure |5] 
we also show projected surface density maps centred around 
the blob position for the 32 x 32 x 128 and 64 x 64 x 256 
runs at t = 3 txH • 

We quantify the time evolution of the blob in Fig- 
ure |6] Here we show the remaining blob mass fraction as 
a function of time, where the material associated with the 
blob is selected such as to satisfy: T < 0.9Tmedium and 
p > 0.64pbiob, as has been done in the previous studies. 
We compare GADGET and moving-mesh AREPO simula- 
tions, performed at three different resolutions (as indicated 
on the legend). While for our highest resolution runs there 
is a broad agreement between GADGET and AREPO for 
t < tKH, blob mass fractions are systematically different 
afterwards. This is exactly when the transition in the mass 
loss rate from the ablation-dominated to the fluid instabil- 
ity dominated regime occurs. In the latter regime, AREPO 
clearly delivers a physically more trustworthy solution. The 
moving-mesh AREPO simulation agrees qualitatively very 
well with the outcome of Eulerian grid codes studied in 
[Agertz et al.| ( [2007| ). However, as noted by [Springel[ ( [20TT1 ), 
there is a small but systematic difference with the moving- 
mesh code delivering slightly higher remaining blob mass 
for t > 1.5 X tKH- As expected, the GADGET simulation 
results agree well with the flndings of Agertz et al. (2007) 



for the other SPH codes. This indicates that the inaccura- 
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Figure 4. Projected surface density maps in units of [M0kpc~^] at times t = ^kh (top row), t = 2tKH (middle row) and t = S^kh 
(bottom row) for GADGET (left-hand panels; a = 1.0, A^ngb = 64) and AREPO (right-hand panels; moving mesh). The thickness of 
the slices is Ay = 100 kpc, and they are centred on yc = 1000 kpc. Dynamical fluid instabilities are responsible for blob shredding on 
the characteristic Kelvin-Helmholtz timescale in the case of AREPO, but they are artiflcially suppressed in the run with GADGET, 
prolonging the survival time of the blob. 



cies we find for GADGET are inherent to the standard SPH 
method, and have prompted a number of works suggesting 
possible improvements to the SPH metho d (e.g. |Price|2008| 



Wadsley et al.|2008||H^ fc Springel|2010| [Saitoh Makino| 
2012[ see also Section T2.1.2| ). 



3.3 



Shocks and fluid instabilities in isolated halo 
models 



3.3.1 Isolated haloes in hydrostatic equilibrium 

To explore more directly the impact of shocks and fluid in- 
stabilities in the context of structure formation, we have de- 
vised a number of idealized test problems involving isolated 
halo models. Here we briefly describe how we set up and ver- 
ify hydrostatic equilibrium configuration of the initial con- 
ditions, which form the backbone for a series of numerical 
experiments discussed in Sections |3.3| and [3^ 

The initial conditions are constructed by populating 
static or live dark matter potentials with gas particles, whose 



positions are drawn randomly. For the dark matter distribu- 
tion we assume a Hernquist profile ( Hernquist|1990 ) so that 
we can easily generate self-consistent models when we use 
live haloes. The gas density profile traces the dark matter 
at large radii but is slightly softened in the centre, i.e. 



(r) 



(27ra3)(a; + a;o)(a; + l) 



(1) 



where Mvir is the virial mass of the system, a is the Hern- 
quist scale length parameter, x = r/a, and xo = 0.01 is 
the softening scale length parameter for the gaseous halo. 
For the simulations without any net rotation, the initial gas 
velocities are set to zero, while for the simulation with non- 
vanishing angular momentum we assign gas velocities such 
that the halo is characterized by the dimensionless spin pa- 
rameter 



A = 



J\E 



1/2 



CM 



5/2 ■ 



(2) 
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Figure 5. Projected surface density maps in units of [M0kpc~^] at t = SIkh for the low (top row) and intermediate resolution simulation 
(bottom row) with GADGET and AREPO. The thickness of the slices is Ay = 100 kpc and they are centred on yc = 1000 kpc. Even at 
our lowest resolution, AREPO captures the dynamical evolution of the cold blob much more accurately. 
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entropy between GADGET and AREPO are within a few per- 
cent throughout the whole halo at the final time in the case 
of static dark matter haloes (see Figure Al in Appendix [A] 
for live haloes see Section 3.5). 



We also observe that even though the initial gas ve- 
locities are zero (in the case with A = 0) some small gas 
velocities develop over time (crgas,3D 30kms~^). This is 
primarily caused by the initial Poisson sampling of gas po- 
sitions, which implies an initial state that is not perfectly 
relaxed. While this numerical artifact could be avoided by 
explicitly relaxing the initial conditions, we note that these 
residual gas velocities do not have any significant bearing 
on our results: the total gas kinetic energy ^kin is less than 
0.005 of either the total potential or the total internal gas en- 
ergy at the final time, as shown in Figure [A2] in Appendix [A] 
It is, however, interesting to note that while the total a gas is 
very similar between GADGET and AREPO, there are some 
systematic differences in the radial profiles of cTgas, which 
are caused by dissipation of gas motions on different spatial 



scales (see bottom panel of Figure Al ). 



Figure 6. The remaining blob mass fraction as a function of time 
in units of ^kh- GADGET and AREPO results at three different 
resolutions are shown, as indicated in the legend. 



where J represents the angular momentum, E is the total 
energy of the halo, and we assume solid body rotation. 

For validation purposes we evolve isolated haloes with 
GADGET and AREPO for 2.45 Gyr with gas self- gravity and 
no radiative losses. These test runs confirm that the gas is 
in very good hydrostatic equilibrium within the dark mat- 
ter potential. The differences in gas density, temperature and 



3.3.2 Radial collapse of cold gas in a static dark matter 
halo 

We now analyze how differences in shock capturing between 
GADGET and AREPO affect the radial infall of gas in dark 
matter haloes, a problem of direct cosmological interest. 
For this purpose we intentionally adopt a set-up as sim- 
ple as possible in order to isolate differences between the 
codes driven by the shock treatment only. We initially set 
up gas in hydrostatic equilibrium and at rest within a static 
Hernquist potential with mass Mvir = lO^^M©, scale length 
a = 176 kpc, and a gas fraction of /gas = 0.17. We consider 
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Figure 7. Radial profiles of gas density, temperature and entropy 
at time t = 0.1 Gyr for GADGET (blue lines) and for AREPO 
(red lines). For each code we show three different resolution runs: 
ATgas = 10^ and r^oft = 14.0 kpc (dotted lines), A^gas = 10^ and 
Tgoft = 6.5 kpc (dashed lines) and A^gas = 10^ and rgoft = 3.0 kpc 
(continuous lines). Vertical black lines with the same style indi- 
cate the softening scales of the dark matter potential (note that 
this does not correspond strictly to the spatial resolution limit 
because the gas is not self-gravitating). The black vertical dot- 
dashed line denotes the virial radius of the system. 
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Figure 8. Radial profiles of gas entropy at time t = 0.1 Gyr for 
GADGET (blue lines) and for AREPO (red hues). For each code 
we show three different resolution runs: ATgas = lO'^ (dotted lines), 
A^gas = 10^ (dashed lines) and A^gas = 10^ (continuous lines), 
which are largely overlapping. These runs have been computed 
by adopting an analytic Hernquist dark matter halo without any 
softening. The black vertical dotted line denotes the virial radius 
of the system. 



both an analytic gravitational potential and a potential with 
a centrally softened cor^ We introduce a small modifica- 
tion in the codes such that gas self-gravity is switched off - 
the gas only feels the static dark matter potential and hydro- 
dynamical forces in a purely non-radiative regime. We then 
artificially reduce the internal energies of gas particles/cells 
so as to displace the gas from the equilibrium solution. The 
newly assigned gas temperature is 4.7 x 10^ K, equal for all 
resolution elements. This test is hence analogous to the well- 
studied Evrard collapse ( Evrard|1988 ) , but it is even simpler 
in nature because we intentionally neglect gas self-gravity. 

As a consequence of the dramatic reduction in its tem- 
perature, the gas will suddenly lose pressure support and 
radially free-fall towards the centre. As the gas collapses, a 
radial shock develops in the centre, steepening while it prop- 
agates outwards and ploughing trough the remainder of the 
outer material which is still falling in. The Mach number of 
the shock varies over the range ^3 — 8, well matched to 
the ID shock tube test problem described in Section [3.1.1| 
Finally, as the shock propagates beyond the virial radius 
of the halo, the gas will reach a new hydrostatic equilib- 
rium solution within the static dark matter potential. In 
Figure [7| we show radial profiles of gas density, temperature 
and entropy computed at time t = 0.1 Gyr after the start 
of the gas collapse for dark matter haloes with softened po- 
tentials. For each code (GADGET, A^gb = 64, a = 1.0: blue 
lines; AREPO: red lines) we perform three runs at different 
resolutions: A^gas = 10^ and rsoft = 14.0 kpc (dotted lines), 
Ngas = 10^ and rsoft = 6.5 kpc (dashed lines) and A^gas = 10^ 
and rsoft = 3.0 kpc (continuous lines). Note that the rsoft 
values indicate the spatial scale over which we smooth the 
analytic dark matter potential. They do not necessarily rep- 



^ We modify the analytic Hernquist potential by convolving it 
with the spline softened potential in the centre, as we describe in 
detail in Section 13. 51 
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resent the minimum spatial resolution of these simulations, 
given that the gas is not self-gravitating and that the dark 
matter halo is rigid. The smoothing lengths of gas particles 
in the central region are of the order of rhsmi ~ 1 kpc in GAD- 
GET for the lowest resolution run with A^gas = 10^, while the 
typical central cell sizes in AREPO are a few hundred pc. 

From Figure |7| it can be seen that there are systematic 
differences in the gas properties predicted by the simulations 
with GADGET and AREPO. In particular, the gas entropy 
distribution is broader in GADGET both in pre- and post- 
shock gas in all three simulations, while for the two lower 
resolution runs there is a slight mismatch in the exact po- 
sition of the shock front and in its strength (similar to the 
findings for the Evrard collapse in Springel|2010a| ) , which is 
minimized in the case of iVgas = 10 particles. Note that the 
differences in shock front position between different resolu- 
tion simulations for a given code are not driven by resolu- 
tion effects but by different spatial softening of the central 
potential, which affects the gas collapse in the innermost re- 
gions. We have checked this explicitly by performing runs 
with A^gas = 10^, 10^ and 10^ for both codes, but this time 
simulating cold gas collapse within an analytic Hernquist 
potential, as shown in Figure [8] In this case the shock prop- 
erties for different resolution runs are almost identical for a 
given code, indicating that the simulation with A^gas = 10^ 
resolution elements is in principle sufficient for capturing 
the shock position accurately. Nonetheless, regardless of the 
resolution, differences between GADGET and AREPO in the 
pre- and post-shock gas persist also in the case of analytic 
Hernquist potentials. 

The differences in entropy content of pre- and post- 
shock gas in GADGET and AREPO are in part due to the 
large shock broadening in SPH, as discussed in Sections |3.1.1| 
and 3.1.2 In fact, if we adopt A^ngb = 32 (which is considered 
the minimum number still permissable in three-dimensional 
simulations) instead of A^ngb = 64, the gas entropy profile in 
GADGET becomes less broad, but is still not as sharp as in 
AREPO, indicating that simulations with a larger number of 
particles are needed in GADGET than in AREPO to recover 
shock features with the same accuracy. 

Moreover, there is another numerical effect leading to 
spatially different entropy generation in GADGET: in the 
converging subsonic part of the flow, artificial viscosity 
(as implemented in GADGET) leads to artificial dissipation 
which increases the entropy in the pre-shock gas. While this 
feature is clearly visible in Figure 40 of |Springel| ( [20T0at for 
the case of the Evrard collapse, here we see that it also en- 
larges the central entropy in GADGET. The reason for this 
is the following: as soon as the gas is brought out of equi- 
librium and starts free-falling towards the centre, the gas 
entropy will be boosted in the central region due to an ac- 
tive artificial viscosity in the converging flow, creating an 
entropy bump that extends up to several kpc (or even sev- 
eral tens of kpc in the case of potentials with large cores) 
away from the centre, even though the shock has not fully 
formed yet at this point. This artificial entropy generation 
is much smaller in AREPO. As the shock forms and propa- 
gates outwards, it will lead to additional physical dissipation 
of much higher magnitude, bringing the entropy profiles of 
GADGET and AREPO into better agreement. Interestingly, 
in the case of the Evrard collapse, the initial difference in 
central entropy profiles is minimized with time due to the 



gas self-gravity (as we checked explicitly by running a sim- 
ulation with exactly the same gas configuration but with 
gas self-gravity and without static dark matter potential), 
while it persists in our test runs even when the new equilib- 
rium solution of the system is reached. The central entropy 
is higher in GADGET by a factor ^ 1.2 and ^ 1.5 for an- 
alytic and softened potentials, respectively. Thus it follows 
that differences in the gas properties due to different spatial 
dissipation of kinetic energy in GADGET and AREPO can 
be aggravated in the case of non self-gravitating gas. 

Note, however, that the difference in the central entropy 
profiles between the two codes goes in the opposite direction 
to what is found in non-radiative simulations of hierarchi- 
cally forming galaxy clusters, where the central entropy is 
higher in mesh-based calculations. This indicates that accre- 
tion shocks during cosmological structure formation do not 
seem to be the likely cause of this central entropy discrep- 
ancy. 



3.3.3 Infall of two gaseous spheres in a static dark matter 
halo 

We now further increase the complexity of the problem by 
considering two gaseous spheres instead of one, collapsing 
into one common static dark matter halo placed in-between 
the two spheres. Each sphere has an initial spatial displace- 
ment from the centre of the halo. This test is similar in 
spirit to a number of previous works which analyzed colli- 
sions of two galaxy clusters in isolat ion (see e.g. [Ricker 
'Sarazin|200T] [Ritchie Thomas] 2 002) [McCarthy et al.|2007| 



Springel Farrar|2007[ [Mitchell et ai.|2009[ [ZuHone|2011[ ), 

but here we devise a cleaner set-up in order to minimize 
additional possible numerical effects (e.g. gravitational N- 
body heating, differences due to gas self-gravity, etc.). As 
in Section [3.3.2| we simulate a static analytic Hernquist 
dark matter halo (with exactly the same parameters) but 
instead of one cold gaseous sphere we generate two identi- 
cal cold spheres separated by 1.2 Mpc along the x-axis, and 
we again neglect any radiative losses and gas self-gravity. 
The gaseous spheres are constructed in the same way as in 
Section [3.3.2| i.e. gas is first set up to be in hydrostatic equi- 
librium within a static Hernquist dark matter potential, and 
then its internal energy is reduced to 4.7 x 10"^ K. For each 
code we perform runs with different gas particle numbers. 



i.e. ATg,, = 10^ 



= 10^ and ATgas = lO"" per sphere. 



Under the gravitational pull from the central dark mat- 
ter potential the two cold spheres collapse towards its centre 
and violently collide. The interesting aspect of this problem 
is that radial symmetry is broken and the gas interaction 
results in much more complicated shock geometries. This is 
illustrated in Figure [9] where we show a time sequence of 
projected gas density maps (first two columns) and mass- 
weighted entropy maps (last two columns) for runs with 
A^gas = 2 X 10^ performed with the two codes. While the 
GADGET simulation shows qualitatively similar gas struc- 
tures, the detailed properties substantially differ. 

Initially, as the spheres start to fall in, the gas is com- 
pressed in the centre of the halo, generating a spherical over- 
density, which is somewhat broader and less peaked in GAD- 
GET, largely due to a poorer effective spatial resolution and 



non- negligible artificial viscosity (see Section 3.3.2). Also, 
during this initial stage more entropy is produced in the 
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Figure 9. Projected surface density maps (first two columns; in units of [M0kpc~^]) and mass-weighted entropy maps (last two columns; 
in internal units) for GADGET and AREPO simulations, showing the collision of two gaseous spheres with A^gas = 2 x 10^ at times 
t = 0.54 Gyr (first row), t = 1.2 Gyr (second row), t = 2.9 Gyr (third row), and t = 15Gyr (fourth row). The plotted spatial domain is 
0.6 X 0.6 X IMpc. 



central region in our SPH calculation. As more gas falls in, 
a shock develops which rapidly assumes a cocoon-like ge- 
ometry elongated perpendicular to the direction of collapse 
(see top panels of Figure |9|. From the gas density maps it 
can be seen that the shock front is narrower and sharper 
in AREPO, whereas in GADGET it has a more splotchy- 
like appearance, caused by kernel averaging. High entropy 
plumes propagating outwards along the ^/-axis are clearly 
visible in the right-hand panels, where the central entropy 
in GADGET within 250 kpc is still higher than in the 
moving mesh simulation. As the cocoon propagates against 
the inf ailing material, gas in the very centre is pushed per- 
pendicular to the X-axis, generating a dense sheet-like region 



(see second row of Figure [9|. Dynamical fluid instabilities at 
the boundary of this dense region induce typical mushroom- 
like morphologies (cap-like in projection). However, even for 
our highest resolution simulations with A^gas = 2 x 10^ par- 
ticles, there are some marked differences between the two 
codes. Mushroom-like features (corresponding to the red- 
orange colours in the density maps) originating at the very 
boundary of the dense central region are more coherent in 
GADGET, while in AREPO they break up and mix more ef- 
ficiently with the surrounding medium. This also leads to 
the more efficient mixing of different entropy gas in the very 
core in the moving mesh simulation, as can be seen from the 
entropy maps. 
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Figure 10. Radial entropy profiles at t ~ 15 Gyr from the start of 
the simulation when the system has reached equilibrium. For each 
code (GADGET: blue lines; AREPO: red lines) three simulations 
with increasing gas particle number are shown: A^gas = 2 x 10^ 
(dotted lines), A^gas = 2 x 10^ (dashed lines), A^gas = 2 x 10^ (con- 
tinuous lines). The vertical dotted line indicates the virial radius 
of the underlying dark matter halo. While for A^gas > 2 x 10^ the 
entropy profiles seem converged for each code, they converge to a 
very different result. 



The differences in the fluid properties at the early stages 
of the simulated system as described above are, however, 
not the only reason why the thermodynamic properties of 
the gas are systematically discrepant between the two codes 
when a new equilibrium state is reached. With time, dense 
shells of gas completely disperse by mixing with the sur- 
rounding material in AREPO, while in the case of GAD- 
GET filaments and blobs of dense gas survive and gradu- 
ally sink back to the centre (see third row of Figure |9|. 
This buoyantly-driven deposition of low entropy material 
in GADGET causes even larger differences between the fi- 
nal entropy distributions, which are illustrated in Figure 
(see also bottom row of Figure [9| . The radial entropy pro- 
files are computed once the system has reached hydrostatic 
equilibrium at time t ^ 15 Gyr from the start of the simula- 
tion. The blue lines denote GADGET results at three differ- 
ent resolutions (dotted lines: A^gas = 2 x 10^, dashed lines: 
ATgas = 2 x 10^, solid lines: A^gas = 2 x 10^), while the red 
lines are for the runs with AREPO. For A^gas > 2 x 10^, 
both codes seem to produce converged entropy profiles, but 
they converge to very different results. While in GADGET 
the entropy profiles steadily decrease towards the centre, 
in the moving mesh code a large entropy core is produced. 
This systematic difference in the central entropy profiles is 
in good agreement with the previous study by Mitche ll et al.| 
( ^2009 J of idealized major merger simulations of haloes in the 
non-radiative regime. 



3.3.4 Generalized blob test: non- radiative case 

We now devise a numerical experiment to capture the evo- 
lution of cold, dense blobs in a more realistic setting, rather 
than a uniform density windtunnel (see Section 3.2.1). For 



this purpose we consider our default isolated halo with a 
static, analytic Hernquist profile for the dark matter com- 
ponent, gas in hydrostatic equilibrium, and we include gas 
self- gravity, but neglect any radiative losses. We additionally 
populate this halo with 10 blobs, with the following prop- 
erties: the gas pressure within blobs is set to be 0.01 of the 
maximum intracluster pressure within 800 kpc radius from 
the centre; the radius of each blob is i^biob = 20 kpc; the 
total mass of all blobs is 20% of the total intracluster gas 
mass, i.e. Mbiobs = 0.2/gasMhaio; the adiabatic index is the 
same for the blobs and for the surrounding gas, 7 = 5/3. 

We place the blobs randomly within a spherical shell 
with cluster-centric distance of 700 kpc and thickness of 
150 kpc, and particles belonging to the blobs are drawn ran- 
domly as well. Apart from positions, we also assign bulk 
velocities to the blobs, while the reminder of the intraclus- 
ter gas is initially at rest. We set all three components of the 
velocity vector for each blob, i.e. Vr (only inward radial ve- 
locity) , V0 , and Vcj) , assuming a random distribution for each 
velocity component starting from a characteristic velocity 
value of 200kms~^. Blob velocity values range then from 
~ 230kms~^ to ~ 510kms~^, with the average velocity of 
all blobs being Vmean ^ 400kms~^. In this way, individual 
blobs will not reach the cluster centre all at the same time, 
giving them more realistic orbits than purely radial ones. Ini- 
tially, the blobs are roughly in pressure equilibrium with the 
surrounding gaseous halo. We perform this numerical exper- 



iment at two different resolutions: A^^ 



10^ 

^5 



lO'' 



(low resolution simulation), and A^gas = 10^, A^biob = 10^ 
(higher resolution simulation). 

In Figure [TT] we show the time evolution of the dense 
blobs moving through the isolated halo, for the higher reso- 
lution run. In the top panels, the projected surface density 
map is plotted at the initial time, where both the intraclus- 
ter gas structure and the blob properties are exactly the 
same in the two codes. Initially, for t < 1 Gyr, the blobs are 
moving on almost identical orbits in GADGET and AREPO, 
and they also have very similar morphologies. As the blobs 
start approaching the inner cluster region, well-defined bow 
shocks develop ahead of each blob and ram-pressure strip- 
ping ablates the blobs. Additionally, the Kelvin-Helmholtz 
and Ray leigh- Taylor instabilities arise which tend to dis- 
rupt the blobs on a characteristic timescale of several Gyrs, 
as discussed in Section 



3.2.1 



(but note that here gas is self- 
gravitating). The lower rate of ram pressure stripping and 
suppression of fluid instabilities in GADGET has a signifi- 
cant effect not only on blob morphologies, but also on their 
orbits. At t = 1.37 Gyr, as illustrated in the middle panels, 
it is still possible to cross-identify each blob in GADGET 
with the respective blob in the AREPO run, but the blobs 
in GADGET have lost less material and thus appear denser 
and some of them are closer to the centre. 

This different loss of blob material leads to systemati- 
cally diverging orbits due to higher buoyancy and dynam- 
ical friction forces acting on each blob in GADGET. This 
can be clearly seen in the bottom panels where the blobs 
in GADGET are markedly more concentrated and have es- 
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Figure 11. Projected surface density maps in units of [MQkpc""^] at the initial time (top panels), dX t = 1.37 Gyr (middle panels), and 
at t = 2.65 Gyr (bottom panels) for GADGET (a = 1.0, A^ngb = 64) and AREPO. The plotted spatial domain is 2 x 2 x 2Mpc, so as to 
encompass initially the whole halo whose virial radius is R200 = 755 kpc. While in the beginning the orbits and the morphologies of the 
blobs are very similar in the two codes, for t > 1 Gyr they begin to diverge, with blobs sinking to the centre faster in GADGET. 
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Figure 12. Spatial distribution of gas at t = 1.37 Gyr that was contained in the blobs at time t = 0. The plotted spatial domain is 
2 X 2 X 2Mpc (for details on map-making see the main text). Most of the material stays confined within blobs in the case of GADGET, 
with a relatively small fraction lost to their wake. The spatial distribution of blob material is markedly different in AREPO, showing 
much more stripped gas that creates prominent tails which extend up to several 100 kpc. 
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Figure 13. Projected mass-weighted vorticity maps in units of [km s~-'^kpc~"'^] (absolute value of the z-component only) for a run with 
GADGET (left-hand panel) and with AREPO (right-hand panel). The maps are computed at time t = 1.37 Gyr from the start of the 
simulation. The thickness of the projected region is Az = IMpc. In the wake of the bow shocks produced by the moving blobs turbulence 
is generated. The spatial extent of the turbulent wakes is significantly larger in AREPO. 



sentially all fallen to the innermost cluster region, while in 
AREPO they are at much larger cluster-centric distances be- 
ing gradually eroded. Moreover, given that the blobs in the 
moving mesh calculation are ablated more efficiently and 
that therefore the dynamical friction exerted on the blobs 
is lower, they have higher velocities when passing at the 



pericentre and thus can reach larger distances after the first 
passage, as visible in the bottom panels of Figure [TT] As 
the gas is self- gravitating, the blobs are also subject to tidal 
stripping when passing close to the innermost regions which 
contributes to the ablation of the blob material. 

To illustrate more clearly the differences in mass loss of 
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in Figure [12] we show 
1.37 Gyr (correspond- 



the blobs in GADGET and AREPO, 
projected maps of gas material at t = 
ing to the middle panels of Figure 1 1 ) which initially be- 
longed to the blobs. To compute these maps in the case of 
GADGET, we show integrated mass along the line-of-sight 
(Az — 2 Mpc) for all particles initially contained in the 
blobs. In the case of AREPO, we use a tracer field to fol- 
low spreading of fluid elements initially within the blobs: 
each cell is characterized by an additional scalar field Tracer 
which at t = is equal to 1 for all blob cells and else- 
where. The tracer field essentially evolves as a dye cast on 
the moving fluid, and at some time t > the Tracer value in- 
dicates the mass fraction of the material which initially was 
in the blobs for each cell. In the right-hand panel of Fig- 
ure [12] we plot the projected density- weighted tracer field. 
The dynamical range is the same in both panels, ranging 
from the maximum of the projected quantity to 10~^ of this 
maximum value, using a logarithmic colour mapping. As an- 
ticipated, cold dense blobs in GADGET lose much less ma- 
terial while in the run with AREPO the lost blob material is 
significantly more diffuse. In the wake of the inf ailing blobs 
prominent tails develop, extending up to several 100 kpc. 
The mass deposition of blob material is clearly much more 
spatially extended in the moving-mesh code and occurs over 
larger cluster-centric distances, allowing fluids with different 
entropies to intermingle and to affect host halo properties 
on larger scales. 

The formation of bow shocks in front of the moving 
blobs also implies that vorticity will be generated due to 
the baroclinic term. To explore this issue in more detail, in 
Figure [13] we show vorticity maps for the run with GAD- 
GET (left-hand panel) and for the simulation with AREPO 
(right-hand panel) at t = 1.37 Gyr (corresponding to the 
middle panels of Figure 11). The vorticity maps are con- 
structed by first evaluating projected mass- weighted veloc- 
ity maps for the x and ^/-components, and then by taking 
the absolute value of the z-component of the curl. Note that 
while the positions and the structure of the bow shocks are 
rather similar in the two codes, as expected, there is a strik- 
ing difference between the spatial extent of the high vor- 
ticity regions produced in the wake of the bow shocks (de- 
noted with green colours), where turbulent motions should 
be generated. For example, focusing on the right-most blob 
centred on x ^ 250 kpc and y ^ —100 kpc in projection, the 
mean projected vorticity value in its wake is a factor of ^ 2 
higher in the moving mesh calculation. This clearly suggests 
that the suppression of vorticity generation in GADGET will 
have an impact on the level of turbulence injected by curved 
shocks that are associated with galaxy or subhalo motions 
through the intracluster medium (ICM), producing a bias 
in the amount of non-thermal pressure support in galaxy 
clusters (see also e.g. Dolag et al. 
lapichino et al.||2011 ). 
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As the blobs orbit at larger cluster-centric distances in 
AREPO for a longer time and mix with the surrounding 
medium more efficiently than is the case for GADGET, they 
lead to higher entropy generation over a wider range of radii. 
This is shown in Figure [l4l where we compute radial entropy 
profiles for all intracluster gas including the blobs, at the fi- 
nal time t 10 Gyr when the system has reached hydrostatic 
equilibrium. Both in the low and high resolution runs with 
AREPO the gas entropy profile is significantly higher in the 
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Figure 14. Radial entropy profiles at time t ~ 10 Gyr after the 
start of the simulation, at which point all the blobs have been dis- 
rupted and the system has reached a new equilibrium. For each 
code (GADGET: blue lines; AREPO: red lines) two resolution 
simulations are shown: A^gas = lO'^ and A/'biob = 10*^ (dashed 
lines), A^gas = 10^ and A/'biob = 10^ (continuous lines). Grey 
curves with the same line-styles represent the initial entropy pro- 
files of these simulations. Additional green line shows AREPO 
simulation with Ngas = 10^ and A/'biob = 10*^ , but with the fixed 
gravitational softening as in GADGET. AREPO runs with adap- 
tive and fixed softenings for the gas produce very similar entropy 
profiles. Vertical lines with the same line-styles denote the gravi- 
tational softenings and vertical dotted line the virial radius of the 
system. The interaction of the moving dense blobs with the intra- 
cluster medium leaves systematically different imprints on the gas 
entropy profiles when simulated with AREPO or with GADGET. 



inner regions up to r ^ 100 kpc. Instead, in GADGET, the 
entropy content of the dense blobs is increased less as they 
sink towards the centre, such that they settle on a lower 
adiabat, corresponding to smaller radii. We also note that 
these systematic differences in entropy profiles are entirely 
due to the different hydro solvers employed by GADGET 
and AREPO, and not due to the different choice of gravita- 
tional softenings for gas (fixed in GADGET and adaptive in 
AREPO, but with a floor set equal to the GADGET value). 
In fact, in Figure [M] the green line shows the entropy profile 
obtained from an identical low resolution AREPO simulation 
where instead of an adaptive a constant gravitational soft- 
enings is used for the gas, in the same way as in GADGET. 

It is also instructive to analyze how the entropy profiles 
of the intracluster medium evolve with time. In Figure [14] we 
show radial entropy profiles at the initial time (grey lines) 
for both resolutions. The initial entropy profile of the run 
with a higher particle number extends further inwards due 
to the better spatial resolution, while some differences in 
the outer regions are due to different positions of the blobs 
which are drawn randomly. The kink in the initial entropy 
profile for r > 500 kpc is caused by the blobs which popu- 
late this region and have low entropy content with respect 
to the surrounding ICM. During the first 0.5 Gyr, the en- 
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tropy profiles in GADGET and AREPO do not ciiange no- 
ticeably, remaining nearly identical to each other and to the 
values prescribed by the initial conditions. At 1 Gyr, the 
gas entropy starts rising in the very centre in the run with 
AREPO (it is higher by a factor of < 3 with respect to the 
initial value) , while it takes almost twice as much time in the 
run with GADGET before the central entropy rises by the 
same amount. At that time, however, the central entropy 
profile in AREPO is already about a factor of 6 higher than 
it was initially, and this systematic difference in the central 
entropy values persists until the final simulated time. 

At time t ^ 2 Gyr another interesting process occurs: 
as several blobs reach the cluster core region for the first 
time, dissipation of their kinetic energy leads to a gradual 
expansion of the whole gaseous halo (compare middle to 



bottom panels of Figure 11) which needs to readjust itself 
to find a new equilibrium solution within the static dark 
matter potential. 



Interestingly, in a recent paper by Vazza et al. (2011) 



non-radiative cosmological simulations performed with the 
GADGET and ENZO codes show that infalling satellites in 
GADGET sink to smaller radii and are characterized by lower 
entropy content than is the case for ENZO, confirming the 
importance of the physical processes discussed here in the 
full cosmological setting. This indicates that at least part 
of the systematic difference in central entropy values found 
between SPH and grid codes in the Santa Barbara cluster 
comparison project ( Frenk et al.|1999 ) is due to the different 
treatment of stripping and mixing of the cold material from 
the infalling satellites. 



3.4 Radiative gas cooling in isolated haloes 

3.4-1 Radiative gas cooling in non-rotating haloes 

We now investigate the time evolution of an isolated 10^^ M© 
mass halo subject to radiative cooling and star formation. 
We first simulate a halo without any net rotation and com- 
pare gas and stellar properties of this system between the 
two codes. The results illustrated here are for haloes with 
static dark matter potentials, but we find very similar results 
if we consider live haloes instead. 

In Figure^] we show star formation rates as a function 
of time for our simulated object consisting of A^gas = 10^ or 
Ngas = 10^ particles/cells, calculated either with GADGET 
or AREPO. Over the whole simulated time-span of 10 Gyr 
the star formation rates are very similar in the two codes. 
This indicates that not only do the implemented gas cooling 
rates match very well, but also that the sub-grid model for 
star formation and supernova feedback is consistent between 
the codes, even at a relatively low numerical resolution (a 
similar conclusion is reached for mergers of isolated galaxies 
by |Hayward et al.|[2Q12J . 

To demonstrate more clearly that gas radiative cool- 
ing and star formation proceed in a very similar manner 
in our isolated haloes in the absence of any net rotation. 



in Figure 16 we plot the total mass of cold baryons (stars 
and gas with entropy < 10^ in internal units) as a function 
of time for GADGET (blue lines) and AREPO (red hues), 
at two different resolutions, i.e. A^part = 10^ (dotted lines) 
and A^part = 10^. For the same number of SPH particles as 
cells used in AREPO the amount of cold baryons matches to 
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Figure 15. Time evolution of the star formation rate of a 
Mvir = lO-*^^ M0 isolated halo which radiatively cools and has 
negligible spin. Illustrated are low (A^gas = 10^; dotted lines) and 
high (A^gas = 10^; continuous lines) resolution simulations with 
GADGET (blue) and AREPO (red). The agreement in star for- 
mation rates is very good over the whole simulated time-span. 




t [Gyr] 

Figure 16. Time evolution of the cold baryonic mass (stars and 
gas with entropy < 10^ in internal units) for a Mvir = 10"*^^ Mq 
isolated halo which radiatively cools and has negligible spin. Illus- 
trated are low (A^gas = 10^; dotted lines) and high (A^gas = 10^; 
continuous lines) resolution simulations with GADGET (blue) 
and AREPO (red). In the inset, the ratio of Mcoid values found 
in AREPO and GADGET is shown for both numerical resolu- 
tions. The agreement of the amount of cold baryons formed is 
extremely good. 
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Figure 17. 2D histogram of the gas entropy as a function of radial 
distance from the centre for GADGET (blue) and AREPO (red) 
at time t = 5 Gyr for A^part = 10^ (dotted hues) and A^part = 10^ 
(continuous lines). The gaseous halo is allowed to radiatively cool 
but there is no net rotation. The vertical dotted lines indicate the 
gravitational softening used in the GADGET run, which corre- 
spond to the floor values of the adaptive softenings in the moving 
mesh calculations. 



within a few percent between the two codes, as is evident 
from the inset in the plot where we show the ratio of Mcoid 
values found with our moving-mesh code and with GAD- 
GET. In the case of live haloes we find that the total mass 
of cold baryons exhibits the same level of agreement. 

Furthermore, in Figure [TT] we show a 2D histogram of 
intracluster gas entropies as a function of cluster-centric dis- 
tance, after 5 Gyr from the start of the simulation for the 
low resolution and high resolution runs. The distribution of 
gas entropies is almost identical in the two codes for r > rsoft 
(note that below rsoft the simulation results are not trust- 
worthy due to the limited gravitational resolution on these 
scales). It can be seen that there is some difference in the 
gas entropy close to the virial radius of the system, which is 
caused by different boundary conditions (vacuum for GAD- 
GET and a uniform low resolution grid for AREPO). 

These results are in line with the expectation that gas 
cooling and condensation in this simulated system is de- 
termined entirely by the gas properties at the cooling ra- 



dius (|Bertschinger||1989| [White Frenk||1991| |Hernquist"fc 



Springel||2003 ) which corresponds very closely between the 
simulation codes. 



3.4-2 Radiative gas cooling in rotating haloes 

Even though gas cooling and star formation proceed in a 
remarkably similar way in GADGET and the moving-mesh 
code for haloes with vanishing spins, this is not guaranteed 
to remain the case once some degree of rotation is included. 
We therefore simulate exactly the same isolated haloes as 
in the previous section, but imposing a certain level of gas 
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Figure 18. Time evolution of the cold baryonic mass (stars 
and gas with entropy < 10^ in internal units) for a Mvir = 
lO"*^^ M0 isolated halo which radiatively cools and has a spin of 
A = 0.4. Illustrated are low (A^gas = lO'^; dotted lines) and high 
(Agas = 10^; continuous lines) resolution simulations with GAD- 
GET (blue) and AREPO (red). In the inset, the ratio of Mcoid 
values found in AREPO and GADGET is shown for both nu- 
merical resolutions. The agreement between the amount of cold 
baryons formed is poorer for A^gas = 10"^, but improves for the 
high resolution run, with Mcoid in AREPO being ~ 10 — 15% 
higher. 



rotation within the static dark matter potential. To highlight 
the effect, we use a large spin parameter equal to A = 0.4. 

Figure ^] shows the time evolution of the total mass 
in cold baryons (stars and gas with entropy < 10^ in inter- 
nal units) for the rotating haloes simulated with GADGET 
(blue curves) and AREPO (red curves) at different resolu- 
tions. Comparing Mcoid with our findings from Figure [16] for 
non-rotating haloes indicates that overall less gas cools from 
the hot phase if the gas spins. This effect is not surprising 
given our initial conditions. Even though the gas density 
and temperature distribution are initially identical, in the 
simulations where there is considerable spin the gas will be 
subject to centrifugal accelerations, preventing it from col- 
lapsing radially. The (partial) centrifugal support will tend 
to reduce the gas densities and hence the cooling rates. 

More importantly, from Figure [18] it can be seen that 
there is poorer agreement in the amount of cold baryons 
between GADGET and AREPO for a rotating gaseous halo. 
The discrepancy is larger for the low resolution run with 
Ngas = 10^, where the final Mcoid value after 10 Gyr is about 
30% higher in the moving-mesh code. The reason for this 
discrepancy is twofold: at low resolution, GADGET under- 
estimates the cooling rate somewhat, while AREPO overes- 
timates it. The SPH result turns out to be more stable at 
poor resolution than the mesh-based calculation. Here it is 
advantageous for SPH that even for few particles a clearly 
defined phase boundary between cold and hot gas is main- 
tained (in fact, a 'pressure blip' in SPH leads to a sampling 
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Figure 19. 2D histograms of gas entropy as a function of radial distance from the centre for GADGET and AREPO dX t = 5 Gyr for 
the lO"*^^ M0 mass halo which radiatively cools and spins at low (A^gas = 10^; left-hand panel) and high (A^gas = 10^; right-hand panel) 
resolution. The vertical dotted line is the gravitational softening used in the GADGET run, which correspond to the floor value of the 
adaptive softenings in the moving mesh calculations. 



gap at this boundary) whereas in AREPO this boundary 
is blurred at low resolution, causing slightly elevated cool- 
ing. This trend is also confirmed by the gas radial veloci- 
ties which are least negative in the low resolution GADGET 
run and most negative in the low resolution AREPO simula- 
tion. The radial velocities systematically differ within a few 
100 kpc and throughout most of the simulated time-span. 
For A^gas = 10^, the radial velocities obtained with the two 
codes are much closer, which also translates into a better 
match of the amount of cold baryons, with Mcoid higher by 
- 10 - 15% in AREPO. 

At least part of the remaining difference in Mcoid be- 
tween GADGET and AREPO is driven by the interaction 
between the gas which is cooling towards the centre and 
the cold gas that is already in the disk. This is shown in 
Figure where we plot 2D histograms of the gas entropy 
as a function of cluster-centric distance at t = 5 Gyr for 
the low (left-hand panel) and high resolution simulations 
(right-hand panel). In the moving-mesh code the gas sur- 
rounding the central disk has a range of entropy values, with 
some fluid elements exhibiting very large entropies due to 
the shock in the immediate vicinity of the disk. Conversely, 
in the runs with GADGET and especially at low resolution, 
there is a clear gap between cold material in the disk and 



hotter gas in the halo (see e.g. Agertz et al.||2007 Springell 



2010b). Similarly, as described in Section 3.3.2| SPH par- 
ticles are affected by artificial viscosity in the converging 
part of the flow which can slightly offset cooling losses. Note 
that for A^gas = 10^ the difference in gas entropy structure 
around the disk is lower between the codes, as expected, 
given that unwanted artificial viscosity effects are reduced 
and that the gap between the cold and hot gas phase due to 
repulsive pressure forces is smaller. 

From Figure [19] is it also evident that the extent of the 



cold disk is different between the codes, especially at low 
resolutions. To quantify this important effect, in Table [l] 
we summarize the main properties of the forming disk. For 
N^^s = 10^ the half- mass radius of the gaseous disk in GAD- 
GET is almost a factor of two smaller than in the higher 
resolution run, while in the moving mesh code i?gas,HM is 
70% of the value we obtain with A^gas = 10^. This indi- 
cates that the convergence rate of the gas disk size is slower 
in the case of GADGET, due to spurious transfer of angu- 
lar momentum from the cold to the hot phase ( Okamoto] 



et al. [2003 ) and due to the artificial viscosity in the case of 



poorly resolved disks (note, however, that the total angu- 
lar momentum is manifestly conserved in GADGET). In the 
case of the stellar disks they are essentially not resolved in 
our low resolution runs (i^stars,HM ^ Tsoft = 14 kpc), while 
the half-mass radius is ~ 40% higher in AREPO at higher 
resolution. Contrary to GADGET, total angular momentum 
conservation is not automatically guaranteed in the moving 
mesh code, particularly for disks resolved with only a small 
number of cells. Nonetheless, it is reassuring that i?gas,HM 
increases with higher resolution in AREPO (attesting that 
spurious angular momentum transport inwards is probably 
small) and that the value of i?gas,HM obtained with both 
codes for A^gas = 10^ is relatively close. This indicates that 
the differences in the disk sizes in this numerical experi- 
ment are at least partly driven by resolution effects, and in 
particular by the slow convergence rate of GADGET sim- 
ulations. Comparing the resolution requirements to obtain 
good agreement between the codes here with Section |3.4| 
we note that in the case of cooling gaseous haloes with a 
significant spin more resolution elements are needed to fol- 
low accurately the thermo-dynamical interaction between 
the central cold disk and the inflowing hotter gas. 
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code 




-Rgas,max 

[kpc] 


^gas,HM 

[kpc] 


^gas,HM 

[lOiOMo] 


^stars,HM 

[kpc] 


^stars,HM 

[lOiOMo] 


Mcold 

[lOio Mo] 


GADGET 




40.0 


30.3 


27.0 


13.5 


215.3 


484.6 


GADGET 


10^ 


79.0 


54.7 


9.8 


6.8 


238.3 


496.2 


AREPO 




79.0 


47.4 


27.7 


18.1 


267.9 


591.2 


AREPO 


10^ 


85.0 


67.1 


12.8 


9.7 


263.1 


551.8 



Table 1. Gaseous and stellar properties of the disk at t = 5 Gyr which forms in an isolated IO^^Mq halo which radiatively cools and 
rotates. For simulations with two different resolutions with GADGET and AREPO, we list the maximum radial extent of the gaseous 
disk, its half-mass radius and the total mass within the half-mass radius, in the third, fourth and fifth columns, respectively. Stellar 
half-mass radius and the total stellar mass enclosed within it are shown in the sixth and seventh columns. In the eighth column the total 
mass of cold baryons is given, being equal to 2 x (Mgas,HM + ^stars,HM)- All gas particles/cells with entropy less than 10^ (in internal 
units) are assumed to be part of the disk, while for the stellar disk we consider all star particles that form in the simulated volume. 



code A^gas A/'blob ^gas,max ^gas,HM ^gas,HM ^stars,HM ^stars,HM ^cold 

[kpc] [kpc] [10^0 Mo] [kpc] [lO^^Mo] [10^0 Mo] 



GADGET 


10^ 


102 


22.0 


16.3 


12.9 


13.5 


184.4 


394.8 


GADGET 


10^ 


10^ 


45.0 


17.9 


7.5 


15.8 


221.1 


457.1 


GADGET_EOS 


10^ 


103 


24.0 


7.6 


154.4 


/ 


/ 


308.9 


GADGET 


10^ 


10^ 


56.0 


32.9 


5.4 


26.6 


235.1 


481.0 


AREPO 


10^ 


102 


76.0 


37.8 


16.8 


17.8 


199.9 


433.5 


AREPO 


10^ 


10^ 


71.0 


43.0 


17.3 


21.8 


261.3 


557.2 


AREPO_EOS 


10^ 


103 


76.0 


13.4 


244.7 


/ 


/ 


489.5 


AREPO 


10^ 


10^ 


75.0 


41.5 


12.7 


32.7 


273.2 


571.9 



Table 2. Gaseous and stellar properties of the disk at t = 5 Gyr that formed in an isolated 10-*^^ Mo halo which radiatively cools, rotates 
and contains 10 orbiting substructures. For simulations with three different resolutions with GADGET and AREPO we list the maximum 
radial extent of the gaseous disk, its half-mass radius and the total mass within the half-mass radius, in the fourth, fifth and sixth columns, 
respectively. Stellar half-mass radius and the total stellar mass enclosed within it are shown in the seventh and eighth columns. In the 
ninth column the total mass of cold baryons is given, being equal to 2 x (Mgas,HM + -^stars,HM)- All gas particles/cells with entropy 
less than 10^ (in internal units) and within i^gas,max are assumed to be part of the disk, while for the stellar disk we consider all star 
particles that form in the simulated volume (including the blobs). We also list the gaseous disk properties in two additional simulations, 
denoted by EOS, performed at intermediate resolution. In these simulations we adopt the standard sub-grid model for star formation, 
with the dense, cold gas lying on the effective equation-of-state, but we prevent any spawning of star particles. 



3.4-3 Generalized blob test: radiative case 

Here we consider the evolution of cold, dense blobs embed- 
ded in a galaxy cluster as in Section |3.3.4| but now the 
whole system is allowed to radiatively cool. To add an addi- 
tional layer of realism the blobs are constructed to be similar 
to cosmological substructures: they are equipped with their 
own dark matter halo, and stars may form out of their gas 
during the simulated time-span. More specifically, each blob 
is represented by a live Hernquist dark matter halo and gas 
in hydrostatic equilibrium. The virial mass of each blob is 
2 X lO^^Mo, the scale length parameter is a = 41.4 kpc, and 
the gas fraction amounts to /gas = 0.17. As before, we popu- 
late our default lO^'^Mo halo (that has a static dark matter 
potential) with 10 identical substructures. The procedure for 
assigning blob positions and velocities is the same as in Sec- 
tion [3]3]4] but here we use a different random number seed 
which leads to the different initial positions and velocities 
with respect to Section |3.3.4| The positions of blob centres 
are in the range of ^ 650 to ^ 750 kpc, while the character- 
istic blob velocities range from 200 to 500kms~^. Also, in 
Section |3.3.4| the gas in the halo is initially at rest while 
here the gas has considerable angular momentum, which 
contributes to the relative velocity between the blobs and 



the gas. We simulate the evolution of this system at three 
different resolutions: A^gas = 10^, A^biob = 10^ (low resolu- 
tion), A^gas = 10^, A^biob = 10^ (intermediate resolution), 
and A^gas = 10^, A^biob = 10^ (high resolution), with gas 
self-gravity, and with cooling and the sub-grid model for 
star formation. 

Additionally, we perform two simulations at intermedi- 
ate resolution where the standard sub-grid model for star 
formation is included as well, but where we simply pre- 
vent any star particles from being spawned out of the star- 
forming phase, denoted by EOS. This simulation set-up is 
particularly useful for following the thermodynamical evo- 
lution of cold and hot gas for many Gyrs without the dense 
cold gas being subject to fragmentation that would likely 
occur in pure cooling runs. For numerical experiments with 
AREPO, as a default choice, we have not used mesh refine- 
ment and de-refinement. However, we have performed ex- 
tra runs at intermediate resolution where we adopt a de- 
/refinement strategy to limit the mass range of gas cells 
(within a factor of 2 of the gas particle mass in the matching 
GADGET run) which automatically imposes a narrow range 
of stellar masses as well. For these runs we have increased 
the number of gas cells/particles in each blob to 2.5 x 10^, 
such that the gas particle/cell mass in blobs is exactly iden- 
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tical to the one in the parent halo (which is optimal for 
our de-/refinement method). These test runs recover very 
closely all of our results with the default set-up, confirming 
that N-body heating effects (e.g. due to the spectrum of star 
particle masses) are not very important here. 

In Figure [20] we show projected surface density maps at 
t = 0.9 Gyr, t = 2.3 Gyr, and t = 6 Gyr for simulations where 
we use the sub-grid model for star formation, but prevent the 
spawning of star particles. Initially the evolution of the blobs 
proceeds in a very similar fashion in simulations with GAD- 
GET and AREPO. However, already after less than a Gyr 



(see top panels of Figure 20 ) blobs in the moving mesh code 
are much more affected by ram pressure and dynamical fluid 
instabilities, which cause efficient gas stripping. As the blobs 
reach the very inner regions, they interact with the forming 
disk. In GADGET simulations, the blobs have a significantly 
more damaging effect on the disk, simply because they are 
less gas depleted. This is clearly visible in the central pan- 
els of Figure |20| where in GADGET feeble spiral features 
are present, while in AREPO the disk is more extended and 
spiral arms are well developed. Additionally, in the moving- 
mesh simulation a large fraction of the stripped material is 
deposited in the main halo, and is gradually accreted onto 
the central disk, promoting its growth. Blobs in GADGET 
are more coherent and eventually lose their angular momen- 
tum due to hydro-/dynamical friction. As they merge with 
the central disk, an ellipsoidal structure is formed (see bot- 



tom panels of Figure 20). Note that due to the static dark 



matter halo the dynamical friction force arises only due to 
the gas in the main halo. Thus, we expect the blobs will 
lose their angular momentum on an even faster time-scale 
in simulations with live dark matter haloes, as is the case 
for cosmological runs. After 6 Gyrs, the difference between 
GADGET and AREPO simulations is significant: while in the 
latter case an extended disk forms, with gaseous spiral arms 
reaching up to 60 kpc away from the centre, in GADGET we 
are left with a flattened, amorphous blob. 

In numerical experiments where we do not suppress star 
particles from forming, we find very similar results. In Fig- 
ure[2l] we show projected stellar density maps at t = 6 Gyr 
using our intermediate resolution simulations (which both 
in resolution and in time match the bottom panels of Fig- 
ure 20). Clearly, the stellar distribution in GADGET is more 



centrally concentrated, forming a small disk. On the con- 
trary, a well-defined extended stellar disk assembles in the 
moving-mesh calculation, endowed with a central bar. We 
quantify the properties of gaseous and stellar disk in Table|2] 
The convergence rate of the gas disk size in GADGET is even 
slower than in the simulations without blobs, because addi- 
tionally to the spurious transfer of the angular momentum 
the damaging effect of the blobs is the highest in the low res- 
olution runs, where gas stripping is largely suppressed. This 
leads to larger systematic difference in disk sizes between 
the two codes with resolution, but for the highest resolution 
runs AREPO disks are only about ^ 30% larger. 

The interaction of cold blobs with the surrounding 
medium influences the global star formation rates as well. In 
Figure [22l we compare the star formation histories in simula- 
tions without blobs (left-hand panel; see Section 3.4.2 ) with 



cially for the low resolution run due to the more efficient 
gas cooling and thus larger amounts of cold gas available 
for star formation. For the high resolution run the differ- 
ences in star formation rates between the codes are smaller, 
given that difference in the amount of cold gas between the 
codes amount to ^ 30%. The difference in the star forma- 
tion rates between the codes becomes more pronounced once 
we include the blobs. Initially, as can be seen from the inset 
plot where we show the star formation rate from the blobs 
and from the disk separately, the star formation rate within 
blobs contributes more than 50% to the total star forma- 
tion rate, and it is very similar in both codes. This induces 
a much higher peak in the total star formation rate in the 
right-hand panel with respect to the simulations without 
blobs. For t > 1 Gyr the star formation rate in the AREPO 
blobs dramatically drops, and at ~ 2 Gyr it is truncated al- 
together. This is due to ram-pressure stripping and dynam- 
ical fluid instabilities that efficiently remove the gas from 
the blobs. In fact, we can crudely estimate a characteris- 
tic Kelvin-Helmholtz timescale at the initial time using the 



equation from Section 3.2.1 (but note that here gas is self- 
gravitating) , where for the blob density we take the average 
density within the blob radius i?biob (which we varied from 
the scale length parameter to the virial radius) and for the 
surrounding medium density we take the typical halo den- 
sity at the position of the blobs. With these assumptions 
we obtain typical values for txH in the range of 1 — 3 Gyr, 
depending on the blob positions, relative velocities and the 
choice of i?biob- These txH values are comparable to the 
timescale on which star formation within the blobs is sup- 
pressed in AREPO. In GADGET simulations, star formation 
proceeds in the blobs even until 9 Gyrs, albeit at a progres- 
sively reduced rate. For t > 1 Gyr the total star formation 
rate comes mainly from the central regions. The central disks 
that form in the AREPO simulations have higher star for- 
mation rates over many Gyrs, which are typically larger by 
a factor of 2 than the GADGET results. Even in our highest 
resolution moving mesh simulation the star forming disk has 
roughly twice as large amount of cold gas, a difference which 
originates from the interaction of hotter infalling gas with 
the gas in the disk, as described in Section |3.4.2| However, 
also the star-forming gas is distributed over a larger area. In 
fact, while within the half- mass radius of the GADGET disk 
the star formation rate is ^ 40% higher in the moving mesh 
run, outside of it is a factor ^ 2.3 higher, contributing about 
half to the total AREPO star formation rate. The reason why 
cold, star forming gas in GADGET is filling a smaller area in 
the disk and is more confined to the dense arm segments and 
blobs (the so called "string-of-pears" , which can be also seen 



in the middle panel of Figure 20 ), is due to the SPH surface 



the numerical tests where haloes are populated with 10 sub- 
structures (right-hand panel). In the case without blobs, the 
star formation rates are somewhat higher in AREPO espe- 



tension originating at the interface between cold and hot me- 
dia in relative motion. In the right-hand panel of Figure [22] 
we also show an identical intermediate resolution AREPO 
run, but with fixed gas gravitational softening, the same as 
in the matching GADGET run. It can be seen that the choice 
of gas gravitational softening in AREPO does not affect our 
results in any significant way: both the evolution of the SFR 
in the central disk and in the orbiting blobs remains very 
similar, indicating that the differences that we see between 
the two codes are indeed entirely driven by hydrodynamical 
effects. 

These findings have immediate consequences for more 
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Figure 20. Projected surface density maps in units of [M0kpc~^] at times t = 0.9 Gyr (top panels), t = 2.3 Gyr (central panels), and 
t = 6 Gyr (bottom panels) for a Mvir = lO-*^^ Mq isolated halo which rotates, radiatively cools, and has 10 orbiting substructures. The 
thickness of the projection is Az = 2 Mpc. Although we have used our sub-grid model for star formation in these simulations, spawning 
of star particles has been intentionally prevented here. Gas stripping from the orbiting blobs is found to be very different in the two 
numerical techniques. In GADGET, blobs largely survive, and when they interact with the central disk they tend to disrupt it, whereas 
in AREPO, the blobs are more easily shredded and have a less damaging effect on the forming disk. In fact, some of the stripped blob 
material ends up contributing to the extended gaseous disk. 
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Figure 21. Projected stellar density maps in units of [M0kpc~^] at t = 6 Gyr for a Mvir = 10^^ Mq isolated halo which rotates, 
radiatively cools, forms stars, and has 10 orbiting substructures. The maps have been constructed from our intermediate resolution runs 
and the thickness of the projection is Az = 2 Mpc. While in GADGET a centrally concentrated disk forms that is surrounded by a feeble 
thick structure, in the simulation with AREPO, a well-defined extended stellar disk assembles with a central bar. 




t [Gyr] t [Gyr] 

Figure 22. Time evolution of the star formation rate of a Mvir = lO-*^^ Mq isolated halo which radiatively cools and rotates. In the 
left-hand panel we show simulations without blobs, while in the right-hand panel the haloes contain 10 orbiting substructures. In the 
right-hand panel, the curves give results from simulations with three different resolutions, with GADGET and AREPO, using A^gas = 10^ 
(LR), A^gas = 10^ (MR), and A^gas = 10^ (HR) particles/cells. In the left-hand panel, only low and high resolution simulations are shown. 
In the simulations without blobs, the star formation rates are somewhat higher in AREPO, but the relative difference increases when we 
include the blobs. Higher star formation rates in AREPO emanate from the central disk which contains a greater amount of cold, star- 
forming gas which fills a larger area. In the inset plot, we compute the total star formation rate occurring within the blobs (dot-dashed 
and continuous lines) and the one coming from the disk (triple dot-dashed and continuous lines) for our intermediate and high resolution 
simulations with GADGET (blue curves) and AREPO (red curves). In the right-hand panel we also show intermediate resolution AREPO 
run with the fixed gravitational softening for the gas (green lines) which produces very similar SFRs (both in the disk and in the blobs) 
as the simulation with the adaptive softening. 
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realistic astrophysical situations. For example, [Puchweiril 
([20T0I) 



et al. ( |2010 | have simulated a high-resolution sample of 
galaxy clusters with GADGET finding that up to 30% of intr- 
acluster stars form in dense cold blobs - remnants of infalling 
satellites. Suppressed dynamical instabilities in GADGET 
enhance the probability of survival for these dense blobs, 
which can then serve as sites of star formation, thus possi- 
bly over-predicting the number of intracluster stars. In fact, 
if we compute the total mass of cold baryons Mcoid (stars 
plus gas above the density threshold for star formation) in 
two cosmological simulations (for further details see Paper I) 
where in one we prevent the spawning of stars and in the 
other we allow it (the simulations are otherwise identical), 
Mcoid is found to be very similar in GADGET regardless 
of whether stars are formed or not. Instead, in the simu- 
lations with AREPO, we find that Mcoid is systematically 
reduced at lower redshifts if the spawning of stars is pre- 
vented. This demonstrates explicitly that very low entropy 
material formed in GADGET cannot easily be shredded and 
mixed with higher entropy gas (at least in the absence of ad- 
ditional feedback processes such as galactic winds or black 
hole heating), so that Mcoid is preserved (a similar conclu- 
sion has been reached independently by |HeB Springel] 
2011). In contrast, in the simulations with AREPO if star 



formation is switched-off some of the cold material can be 
stripped out of galaxies due to dynamical instabilities, re- 
turning it to diffuse form and lowering Mcoid- Importantly, 
this also implies that the number of stars formed in AREPO 
will be much more sensitive to the characteristic timescale 
for star formation and to the numerical resolution (needed 
to resolve fluid instabilities) than is the case for GADGET. 
In a recent paper by |Agertz et al. (2011 ) it has been shown 
that lower star formation efficiency leads to the production 
of larger disks in cosmological simulations. While this result 
is in line with our findings, it is important to recognize that 
the actual cause is quite different: in the study by |Agertz| 
(2011), differences in the physical modelling of star 



et al. 



formation in cosmological simulations affect the disk sizes, 
while here the cause lies in different accuracies of the hydro 
solvers involved. 



3.5 Gravitational N-body heating 

We now investigate potential artificial gas heating due to 
the Poisson noise induced by the finite number of particles 
in dark matter haloes. [Steinmetz White] ( | 1997 ) outlined 
the analytic theory of this effect, and conffimed it with non- 
radiative and radiative numerical experiments that quan- 
tified gravitational N-body noise present in structure for- 
mation simulations. For an equilibrium system they defined 
a characteristic N-body heating timescale which is propor- 
tional to the cube power of the dark matter velocity disper- 
sion and inversely proportional to the dark matter particle 
mass and dark matter density. Due to this inverse propor- 
tionality to the dark matter density, N-body heating is ex- 
pected to be strongest in the innermost regions of haloes. 
From their analysis it follows that the dark matter parti- 
cle mass adopted in numerical simulations should be lower 
than a critical mass which is of order of few times IO^Mq 
for galaxy clusters (see their equation [10]), otherwise the 
radiative cooling losses may be overwhelmed by spurious 
heating. 
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Figure 23. Radial profiles of gas density, temperature and en- 
tropy at time t ~ lOGyr for GADGET (blue lines) and for 
AREPO (red lines). For each code we show two different reso- 
lution runs: A^gas = 10^ and rgoft = 14.0 kpc (dotted lines), and 
A^gas = 10^ and rsoft = 3.0 kpc (continuous lines). The vertical 
black lines with the same style indicate the softening scale of the 
dark matter potential, whereas vertical dot-dashed lines show the 
virial radius of the system. 
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Figure 24. Comparison of fitting functions for spline-softened 
gravitational accelerations (equation |4] blue dashed lines) with 
actual gravitational accelerations obtained from very high reso- 
lution runs with matching rgoft values (green continuous lines). 
The limiting cases of spline-softened accelerations for small and 
large radii are illustrated with dotted red lines. The bottom panel 
shows the residuals of our fitting functions over the relevant range 
of radii. 



Since nowadays the typical mass resolution of hydro- 
dynamical cosmological simulations has improved dramati- 
cally, with simulated galaxy clusters containing several times 
10^ up to 10^ particles in zoom-in runs ( |Dolag et al. 2009; 



Sijacki et al. 



2009 



Vazza et al. 2010), possible numerical 



artifacts due to a grainy dark matter distribution are rarely 
addressed in the literature (but see e.g. Kay et al. 2000| 
Borgani et al.|2006||Vazza|2011 ). Recently, |Springel| (|2010a| 
has pointed out that mesh based codes are potentially more 
severely affected by gravitational N-body heating than SPH 
codes, due to their better ability to detect weak shocks, 
which in this case is an unwanted feature. 

Here we perform a number of idealized test problems 
aimed at explicitly addressing the N-body heating problem, 
and in particular to understand whether there are system- 
atic differences between GADGET and AREPO in this re- 
spect. Note that intentionally all previous tests (except for 



Section 3.4.3) have been performed either without any dark 



matter component or with a static dark matter potential to 
avoid such possible Poisson-noise imprints. 

We first consider our standard isolated halo, where we 
replace the static dark matter potential with a live dark mat- 
ter halo in which dark matter particles have a characteris- 
tic velocity dispersion for the given mass of the system (but 
there is no net rotation). As in the case of the rigid potential, 
we populate the halo with gas particles in hydrostatic equi- 
librium which are initially at rest. We include self-gravity 
of the gas component and evolve the system non-radiatively 
for 10 Gyr. 

In Figure [23] we show radial profiles of gas density, tem- 
perature and entropy at the final time for both codes at 
two resolutions: A^gas = A^dm = 10^ particles (dotted lines) 



and A^gas = A^dm = 10^ particles (continuous lines). The 
centre of the system is defined by the position of the most 
bound particle. For cluster-centric distances greater than the 
gravitational softening value, rsoft, the entropy profiles of 
GADGET and AREPO agree to within 20% for the low 
resolution runs, and to better than 10% for the high res- 
olution simulation. For r < rsoft, the differences between 
the codes in the entropy profile amount up to a factor of 2, 
with the gas entropy in AREPO being systematically higher. 
An extra test simulation performed with AREPO in its dual 
entropy mode (where entropy instead of energy is explic- 
itly conserved for cells whose Riemann problems with ad- 
jacent cells all have Mach number less than 1.1) resulted 
in identical radial profiles as for the AREPO simulation in 
the standard energy mode. This clearly indicates that the 
small entropy core which develops in AREPO is not caused 
by weak shocks generated by the grainy dark matter poten- 
tial. Rather, central gas particles/cells acquire small scale 
velocities when simulated with live haloes. In the case of the 
moving-mesh code, these velocities lead to fluid mixing and 
temperature equilibration in the innermost regions, while in 
GADGET such mixing does not occur. Note, however, that in 
the gravitating systems, rsoft represents the minimum spa- 
tial scale below which the simulation results are not trust- 
worthy due to discreteness effects. Therefore, provided that 
gravitational softening lengths are chosen conservatively, the 
gas properties of galaxy clusters in an equilibrium configura- 
tion simulated with GADGET and AREPO match very well. 

We now consider a more challenging problem, where we 
simulate the infall of cold gas into dark matter haloes. This 
is the same problem discussed in Section |3.3.2| but here 
we compare outcomes from live versus static dark matter 
haloes. To match numerical experiments with rigid and live 
haloes as closely as possible we adopt the following proce- 
dure. 

For live haloes, the motion of dark matter particles is 
prevented in the code to avoid local deformations of the 
gravitational potential which will not be present in the static 
case. Effectively, in this way the distribution of dark matter 
is "frozen", and live haloes are coarse grained representa- 
tions of static haloes. 

For static haloes, we convolve the analytic Hernquist 
potential with a spline softened potential in the centre, such 
that gas particles feel exactly the same gravitational accel- 
eration as in the case of live haloes, where the central poten- 
tial needs to be softened to minimize effects from two-body 
encounters. For the Hernquist potential, gravitational accel- 
erations are given by 



9{r) - 



' (r + a)2 



(3) 



We modify the gravitational acceleration felt by gas particles 
in the code, viz. 



9{r) 



GMvu 



(4) 



where /lo and hi are free coefficients (to first approximation 
ho ^ hi ^ rsoft)- We determine the value of ho and hi for 
three different resolution runs by fitting equation Q to the 
gravitational acceleration of a live "frozen" halo with 10^ 
dark matter particles simulated three times assuming rsoft 
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Figure 25. Radial profiles of gas density, temperature and entropy at time t = 0.3 Gyr for GADGET (left-hand panels; blue lines: static 
dark matter haloes; green lines: live dark matter haloes) and for AREPO (right-hand panels; red lines: static dark matter haloes; green 
lines: live dark matter haloes). For each code we show three different resolution runs: A^dm = 10^ and rgoft = 14.0 kpc (dotted lines), 
A^DM = 10^ and rgoft = 6.5 kpc (dashed lines) and A^dm = 10^ and rgoft = 3.0 kpc (continuous lines), while A^gas = 10^ is kept fixed. In 
the right-hand panel we also show the entropy profile obtained with the moving-mesh code in the entropy mode for our highest resolution 
simulation (continuous orange line). The vertical black lines with the same style indicate the softening scales of the dark matter potential. 
The black vertical dot-dashed lines denote the virial radius of the system. 
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values appropriate for the numerical resolutions we want to 
investigate. 

In Figure [24] we illustrate the outcome of this method. 
In the upper panel, the green lines denote the radial gravi- 
tational acceleration of the live "frozen" halo with 10^ dark 
matter particles, adopting rsoft = 3, 6.5 and 14kpc (indi- 
cated by numbers 1, 2 and 3), respectively. Blue dashed 
lines are our best fit functions to these simulations using 
equation Q. Dotted red lines indicate the radial acceler- 
ation in the case of the analytic Hernquist profile and of 
the spline-softened potential, which are the limiting cases 
of equation (|4| for large and small radii, respectively. The 
bottom panel of Figure [24] shows how accurate our modified 
gravitational acceleration is as a function of cluster-centric 
distance. 

With the procedure described above we perform simu- 
lations of radial infall of cold gas into live ("frozen") haloes. 
We keep the number of gas resolution elements the same 



and equal to A^gas = 10 (corresponding to rsoft, i 



3 kpc 



for GADGET, while gravitational softenings are adaptive 
in AREPO but with a floor of 3 kpc), while we increase 
the number of dark matter particles from A^dm = lO'^ 
(rsoft, DM = 14 kpc), to 10^ (rsoft, DM = 6.5 kpc) and 10^ 
(^soft,DM = 3 kpc). For each live halo we run a match- 
ing static halo simulation where the number of gas parti- 
cles/cells is kept the same, i.e. A^gas = 10^. In all runs, gas 
self-gravity is neglected and there are no radiative losses. 
In close encounters of gas and dark matter, the effective 
gravitational softening is the maximum between rsoft, gas and 
rsoft, DM, such that rsoft, dm is the relevant length scale of the 
problem. 

In Figure [25] we show radial profiles of gas density, tem- 
perature and entropy at time t = 0.3 Gyr (analogous to 
Figure [?]) . The left-hand panels are for simulations of live 
and static haloes with GADGET, while the right-hand pan- 
els show the results with AREPO. Due to gravitational N- 
body heating there are systematic differences in the central 
gas properties. The central gas density is lower, while the 
central gas temperature and entropy are higher for the live 
haloes, due to spurious transfer of energy from dark matter 
to gas, which heats the gas and makes it expand. In general, 
N-body heating effects are largest for A^dm = 10^ and small- 



est for A^E 



10 , and they are confined to spatial regions 



within rsoft, DM, which is reassuring. We can see by compar- 
ing the left-hand to the right-hand panels that AREPO is 
indeed much more sensitive to spurious gravitational heat- 
ing, as discussed in Springel (2010a), with the central en- 
tropy boosted by two orders of magnitude. In the right-hand 
panel of Figure [25] we also show radial entropy profile ob- 
tained with AREPO in the dual entropy mode (where en- 
tropy instead of energy is explicitly conserved for cells with 
all Mach numbers less than 1.1) with A^dm = 10^ particles 
(orange continuous line). This indicates that at least part 
of the central entropy core is generated by the weak shocks 
which gas cells experience when moving through the grainy 
dark matter potential. 

At the final time t — 2.45 Gyr, when the system has 
reached an equilibrium state, the difference between the 
matching live and static halo runs for r > rsoft, dm is ~ 10%, 
~ 5%, and ^ 5%, respectively, for low, intermediate and 
high resolution in GADGET, while it is somewhat higher in 
AREPO, i.e. - 15 - 20%, 10%, and 5%, respectively. For 



r < rsoft the ratios of central entropy values are < 10, 
< 4 and < 2 in GADGET and < 300, < 150 and < 25 
in AREPO, for low, intermediate and high resolution runs. 
Even though the amount of artificial heating is significantly 
larger in AREPO than it is for GADGET, the departures 
between static and live halo radial profiles always occur 
within rsoft, DM for both codes. Provided that gravitational 
softening values for the dark matter component are chosen 
cautiously, our numerical experiments hence indicate that 
the systematic discrepancy in central entropy values found 
between SPH and mesh-based codes for the Santa Bar- 
bara cluster comparison project ( Frenk et al.|1999] Springel 
2010a) is unlikely to be due to effects from a grainy dark 



matter potential during cosmological structure formation. 



4 DISCUSSION AND CONCLUSIONS 

In this study we have carried out a detailed comparison be- 
tween the SPH code GADGET and the new moving-mesh 
code AREPO on a number of hydrodynamical test problems, 
which are crucial for understanding cosmological simulations 
of galaxy formation. In a purely hydrodynamical regime 
without gas self-gravity or an external gravitational poten- 
tial we have first carried out a set of numerical experiments 
previously considered in the literature, some of which have 
rarely been shown for SPH codes, such as the 2D implosion 
test. We have then focused on idealized non-radiative galaxy 
cluster simulations, specifically aimed towards benchmark- 
ing differences in hydro solvers for problems with shocks and 
fluid instabilities. In simulations where radiative losses were 
included we have analyzed the amount of baryons which cool 
from the hot halo atmospheres in GADGET and AREPO, 
both for rotating and non-rotating haloes. In the former 
case, we have also studied how the central baryonic disks 
form, and how orbiting substructures affect the disk mor- 
phology and the star formation rate. Finally, we have con- 
structed special test problems designed to gauge the effect 
of gravitational N-body heating on the gaseous properties of 
the haloes. Our main conclusions are as follows: 

• While post-shock fluid properties are captured well 
both in GADGET and AREPO in the case of ID shock tube 
tests with high Mach numbers, the shocks are significantly 
broader in GADGET, and substantial post-shock oscillations 
develop, largely because of an inadequate treatment of the 
initial contact discontinuity. AREPO in the moving-mesh 
mode preserves the contact discontinuity much more accu- 
rately, but it is broadened if we employ a static mesh, due 
to larger numerical diffusivity in this case. Note that the 
more accurate treatment of the contact discontinuity in the 
moving mesh code can lead to a larger 'wall heating' effect 
(Rider 2000) than in static grid codes, which tend to wash- 
out at some level the initial start-up errors at the contact 
discontinuity (see also description of the Noh problem in 
Springel[[2010a) . 

• For shocks with complicated geometries in multi- 
dimensions, differences between GADGET and AREPO are 
more striking. Even though the global fluid properties re- 
main similar, the sampling of the fluid properties is much 
noisier in GADGET, and the development of dynamical fluid 
instabilities is inhibited. These problems are not specific to 
GADGET, but are inherent to the standard SPH method as 
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a whole (see also' Agertz et al.|20Q7| Springel|2010 b! and dis- 
cussion in Section 2.1.2). The suppression of fluid instabih- 



ties reduces the amount of entropy generation by mixing and 
artificiafly prolongs the lifetime of gaseous structures which 
are moving through a medium with a different density. Ad- 
ditionally, vorticity generation in the wake of curved shocks 
due to the baroclinic source term is largely suppressed in 
GADGET, which directly impacts angular momentum trans- 
fer by vortices and the level of generated turbulence. 

• These fundamental differences between GADGET and 
AREPO identified in simple hydrodynamical test problems 
affect the properties of gas in more realistic, cosmologically 
motivated simulations as well. Specifically, in non-radiative 
idealized simulations of merging galaxy clusters and in simu- 
lations of isolated haloes with orbiting gaseous substructures 
we find that in AREPO: i) an entropy core is produced in the 
centre due to more efficient fiuid mixing, ii) the gas stripping 
rate from the orbiting substructures is larger, and in) more 
vorticity is produced in the wake of curved shocks. These 
findings are in line with results of previous works which sim- 
ulated similar problems with AMR codes (e.g. [Mitchell et al.| 
2009||Vazza et al.||2011| |Vazza||2011J . Moreover, unphysical 
dissipation of shocks and subsonic turbulence in GADGET, 



as shown by "Bauer & Springel (2012) (see also Section 4.2 



of Paper I) , leads to the heating of the halo outskirts rather 
than of the central regions. 

• In radiative simulations of isolated haloes without any 
net rotation, gas cooling and condensation into stars pro- 
ceeds in a very similar fashion in GADGET and AREPO, 
given that the gas properties at the cooling radius match 
closely. However, for spinning haloes there is a net differ- 
ence in the total amount of cold baryons, which is higher 
in AREPO, and this difference persists at a level of about 
10 — 15% in our highest resolution runs. Baryonic disks 
which form due to dissipative collapse of rotating gas are 
systematically larger in AREPO at low resolution, and the 
convergence rate of the gas disk sizes is higher. 

• In numerical experiments where we follow the interac- 
tion between a forming central disk and 10 orbiting substruc- 
tures in an isolated halo which radiatively cools, the final 
disk morphology is significantly different. While in AREPO 
an extended disk is produced with well developed spiral arms 
and a central bar, in GADGET the disk is more centrally con- 
centrated and amorphous. Orbiting substructures are much 
more efficiently stripped of their gas content in the moving- 
mesh calculation and the material is incorporated into the 
host halo atmosphere. Instead, in GADGET gaseous sub- 
structures are more coherent, thus they lose more angu- 
lar momentum from hydro-/dynamical friction, and when 
passing through the central disk they induce morphological 
transformations . 

• While star formation is more readily truncated in in- 
falling substructures due to gas stripping, extended gaseous 
disks in AREPO have significantly larger star formation rates 
for many Gyrs than is the case for GADGET. 

• Due to its better ability to detect weak shocks, AREPO 
is more sensitive to gravitational N-body heating ( S pringel| 
2010a). This is confirmed by our specifically designed nu- 
merical experiments, which allow us to quantify the magni- 
tude of spurious gas heating for simulations of haloes with 
a finite number of dark matter particles with respect to the 
analytic dark matter potentials. Nonetheless, the spatial ex- 



tent of this artificial heating is reassuringly constrained to 
lie within the gravitational softening length for typical set- 
ups, which is a scale below which simulation results are not 
trustworthy due to resolution effects anyway. 

The numerical experiments presented in this study 
clearly demonstrate that several important shortcomings of 
the SPH solver not only affect idealized test problems but 
are equally detrimental in more realistic setups relevant for 
structure formation, and ultimately in full cosmological sim- 
ulations. This is especially the case because of: i) the compli- 
cated fiows involving multiphase media which are the norm 
in cosmological simulations, and ii) the hierarchical nature 
of structure formation where low mass systems are always 
poorly resolved. In both of these regimes the hydrodynami- 
cal solver of the standard SPH method exhibits the largest 
inaccuracies. On the other hand, our numerical tests confirm 
and significantly extend the findings of Springel (2010a) that 
the moving-mesh code AREPO delivers a physically more ac- 
curate representation of the evolution of inviscid gases. Note 
that in the current work we deliberately kept the modelling 
of the baryon physics at a very simple level, so as to isolate 
the differences between the hydro solvers in as clean a man- 
ner as possible. While inclusion of more realistic feedback 
mechanisms, such as supernovae winds and AGN heating, 
will likely modify the properties of the simulated galaxies 
significantly, it is of prime importance to disentangle nu- 
merical inaccuracies of the hydro solver from the uncertain- 
ties of the feedback physics modelling. It is hence clear that 
cosmological simulations with AREPO have the potential to 
provide a much more realistic description of structure for- 
mation in the Universe (see also Papers I and II) , something 
we will explore in more depth in forthcoming studies. 
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Figure Al. Radial profiles of gas density, temperature and 
ID velocity dispersion at time t = 0.05 Gyr (dotted lines) and 
t = 2.45 Gyr (continuous lines) for GADGET (blue lines) and 
AREPO (red lines) with A^gas = 10^ in a static dark matter halo 
with a Hernquist profile. Vertical black lines indicate the soften- 
ing scales of the gas equal to rgoft = 3.0 kpc. The black vertical 
dot-dashed line denotes the virial radius of the system. 
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Figure A2. Time evolution of the ratio of gas kinetic (dotted 
lines), potential (dashed lines), and internal (continuous lines) to 
the total gas energy in GADGET (blue) and AREPO (red) sim- 
ulations of isolated halos in hydrostatic equilibrium. Gas is self- 
gravitating and represented by A^gas = 10^ resolution elements 
and is evolved within a static Hernquist dark matter halo. 



APPENDIX A: ISOLATED HALOES IN 
HYDROSTATIC EQUILIBRIUM 

Here we show how gaseous spheres placed in hydrostatic 
equilibrium within a static Hernquist dark matter halo 
evolve with time to determine the level of accuracy of our 
initial conditions. The general set-up is described in detail 
in Section |3.3.1| Specifically for the figures presented here 
we assume that the gas is self-gravitating and that it has 
initially zero velocity. 

In Figure [AT] we show gas density, temperature and ID 
velocity dispersion radial profiles at t = 0.05 Gyr (dotted 
lines) and t = 2.45 Gyr (continuous lines) for GADGET (blue 
lines) and AREPO (red lines) with A^gas = 10^ resolution el- 
ements. As discussed in Section |3.3.H due to the Poisson 
sampling of the gas positions in the initial conditions they 
are not perfectly relaxed, which leads to the development 
of small-scale random motions, as evidenced by a non-zero 
gas velocity dispersion. Over time these residual gas veloci- 
ties are dissipated, as can be seen from the bottom panel of 
Figure [Al] This leads to a slight readjustment of the tem- 
perature distribution, which is very similar in both codes. 
Note, however, that the kinetic energy of the random gas 
motions is at all times a very small fraction of both the po- 
tential and the internal energy, as illustrated in Figure |A2] 
amounting to only about 0.2% of the total gas energy after 
2.45 Gyr. 
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